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ADDED  MASSES  AND  FORCES  ON  TWO  BODIES 
APPROACHING  CENTRAL  IMPACT  IN  AN  INVISCID  FLUID 

ABSTRACT 

In  several  papers,  which  will  be  referenced,  a  procedure  based  on  integral  equations 
has  been  developed  and  applied  for  determining  the  interaction  forces  on  two  bodies 
approaching  central  impact  in  an  inviscid  fluid.  The  present  work  was  undertaken  to 
evaluate  the  accuracy  of  the  results  from  that  procedure  by  applying  it  to  a  pair  of  circles 
and  a  pair  of  spheres  with  which  one  could  obtain  solutions,  as  accurate  as  desired,  by  the 
method  of  successive  images.  A  second  purpose  was  to  refine  the  procedure  so  that 
accurate  solutions  could  be  obtained  at  closer  distances  than  heretofore. 

Solutions  by  the  method  of  images,  given  by  Hicks  and  Herman  over  100  years 
ago,  are  not  very  clear,  and  since  we  have  significandy  extended  their  theory  in  the  present 
work,  it  seemed  appropriate  to  include  a  new  derivation  which  we  consider  more  rational. 
The  extensions  of  the  theory  consists  of: 

1)  a  truncation  correction  of  the  infinite  series  of  the  doublet  strengths  for  the  added 
masses  and  their  derivatives,  which  can  then  be  calculated  accurately  with  a  moderate 
number  of  terms  even  when  the  gap  between  the  bodies  is  very  small; 

2)  asymptotic  formulas  for  the  added  masses  and  their  derivatives  at  small  gaps 
which  show  that,  for  circles,  the  derivatives  with  respect  to  a  parameter  asymptotically 
proportional  to  the  square  root  of  the  gap,  are  finite,  and  that  derivatives  with  respect  to  the 
gap  approach  infinity  inversely  as  the  square  root  of  the  gap; 

3)  a  treatment  of  the  case  of  a  circular  cylinder  or  a  sphere,  or  bodies  of  arbitrary 
shape  approaching  a  wall,  showing  that  the  forces  on  the  body  and  wall  are  repulsive  and 
of  equal  magnitude; 

4)  a  treatment  of  uniform  convergence  of  the  series  for  the  added  masses  and  their 
derivatives,  which  shows  that  the  infinite  series  for  the  added-mass  coefficients  converge 
uniformly  for  all  values  of  the  gap  in  the  closed  region  from  zero  to  infinity;  and  that  the 
series  of  the  derivatives  with  respect  to  the  parameter  of  item  2)  converge  uniformly  for  all 
values  of  the  gap  except  zero,  where  the  series  converges  but  its  sum  is  discontinuous. 

In  Part  II,  refinements  of  the  integral-equation  procedure  are  presented.  Two 

problems  required  resolution  for  small  gaps.  One  is  the  sharp  peaks  of  the  kernels  of  the 


i 


integral  equations  and  of  the  source  distributions  on  the  body  surfaces  (the  unknown 
functions  of  the  integral  equations)  in  the  neighborhood  of  the  gap.  The  other  is  that, 
previously,  derivatives  of  the  added  masses,  obtained  by  numerical  differentiation,  were 
inaccurate  and  gave  large  errors  in  the  calculated  forces.  Simple  solutions  of  both  problems 
with  the  latter  still  using  numerical  differentiation,  are  presented.  Applications  to  various 
combinations  of  circle  pairs  and  to  equal  spheres  gave  agreement  with  exact  values  to  six 
decimals  for  the  added-mass  coefficients  and  to  five  for  their  derivatives.  To  illustrate  the 
applicability  of  the  new  procedures  to  a  noncircular  cylinder,  the  two  cases  ot  a  2  to  1 
ellipse  approaching  a  circle  in  the  direction  of  its  major  or  minor  axes  are  also  presented. 

It  is  concluded  that  we  now  have  the  capability  of  obtaining  accurate  results  for 
interaction  forces  for  two-dimensional  forms  and  bodies  of  revolution  in  an  inviscid  fluid. 
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Added  Masses  and  Forces  on  Two  Bodies  Approaching  Central 

Impact  in  an  Inviscid  Fluid 

Introduction 

In  connection  with  the  problem  of  an  ice  mass  approaching  central  im¬ 
pact  with  an  offshore  structure,  it  is  important,  for  design  purposes,  to  esti¬ 
mate  the  magnitude  of  the  interaction  force  on  the  structure  prior  to  and  dur¬ 
ing  impact.  Previous  literature  on  this  problem  was  reviewed  in  papers  by 
Landweber,  Chwang,  and  Guo  [1],  and  more  recently  by  Guo  and  Chwang  [2], 
[3],  the  latter  two  mainly  on  the  broader  problem  of  oblique  impact,  with 
some  contributions  on  central-impact.  In  all  of  these  papers,  the  fluid  was 
assumed  to  be  inviscid  and  the  flow  irrotational,  so  that  the  forces  were  en¬ 
tirely  due  to  the  inertia  of  the  two  bodies  and  the  fluid. 

As  is  shown  in  [1],  for  central  impact,  the  interaction  forces  on  the  bod¬ 
ies  due  to  the  fluid  can  be  expressed  in  terms  of  three  added  masses  and  their 
time  derivatives.  The  added  masses,  at  each  of  a  succession  of  instants,  are 
then  obtained  by  applying  generalizations  of  the  Taylor  added-mass  formula, 
given  by  Landweber  [4],  Landweber  &  Yih  [5],  and  Landweber  and  Chwang  [6] 
to  the  sources  and  doublets  on  or  in  each  body.  The  latter  were  obtained  by 
solving,  numerically  and  simultaneously,  a  pair  of  Fredholm  integral  equa¬ 
tions  of  the  second  kind.  These  were  solved  by  discretizing  each  integral 
equation  into  a  set  of  linear  equations. 

Two  major  difficulties  were  encountered  which  affected  the  accuracy  of 
the  numerical  results  when  the  minimum  distance  between  the  bodies  was 
small.  One  is  due  to  the  property  that  the  integrands  of  the  integral  equations 
peak  sharply  at  points  of  the  body  surfaces  which  are  in  the  neighborhood  of 
the  eventual  contact  point  when  the  bodies  have  nearly  met;  the  other  that 
the  absolute  value  of  the  time  derivatives  of  the  added  masses  approach 
infinity  as  the  bodies  approach  contact.  Consequently,  accurate  results  for  the 
added-mass  derivatives  could  not  be  obtained  by  direct  application  of 
numerical  differentiation. 

Techniques  proposed  for  overcoming  these  two  difficulties,  described 
in  three  papers,  were  developed  for  a  pair  of  circular  cylinders,  because,  for 
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this  case,  results  from  the  integral-equation  solution  can  be  compared  against 
those  from  the  method  of  successive  images,  which  can  be  made  as  accurate 
as  one  desires.  Isaacson  and  Cheung  [7]  resolved  the  peaking  problem  by  in¬ 
creasing  the  number  of  points  in  the  neighborhood  of  the  peak,  but  gave  no 
details  concerning  the  quadrature  formulas  used.  Nor  was  the  problem  of 
numerical  differentiation  of  added  masses  mentioned. 

In  Ref.  [1],  "the  most  accurate  quadrature  formula"  [8],  hereafter  called 
the  MAQF,  which  requires  a  smooth,  cyclic  integrand  and  uniform  intervals, 
was  applied,  with  the  number  of  points  used  increased  as  the  gap  between  the 
bodies  decreased.  This  procedure  was  inefficient,  requiring  a  huge  number  of 
points  at  the  smallest  gaps,  where  the  Cray  supercomputer  had  to  be  used. 
The  added-mass  derivatives  were  obtained  by  using  computer  software  to 
"smooth"  the  added-mass  data  before  the  numerical  differentiation. 

In  Ref.  [2],  the  peak  was  partially  removed  by  a  technique  which  is 
available  when  the  kernel  can  be  integrated  exactly;  see  Miloh  and  Landweber 
[9].  Also,  as  in  [7],  the  points  on  the  circles  were  distributed  so  that  at  least  half 
of  their  total  number  were  concentrated  in  the  small  region  of  the  peaks. 
This  combination  yielded  good  agreement  with  the  results  by  successive 
images,  but  has  the  disadvantage  that  the  procedure  used  to  partially  remove 
the  peak  of  the  kernel  is  not  an  option  for  bodies  of  arbitrary  shape.  To 
overcome  the  second  difficulty,  numerical  differentiation  was  avoided  by  dif¬ 
ferentiating  the  integral  equation  pair  and  solving  the  set  of  four  integral 
equations  simultaneously  for  the  source  strengths  and  their  derivatives.  The 
kernels  of  these  new  integral  equations  and  their  solutions  for  the  source- 
strength  derivatives  have  much  stronger  peaks  as  the  gap  approaches  zero,  so 
that  four  Gauss  quadrature  formulas,  each  of  order  40  were  used  to  obtain 
accurate  results  at  small  gaps.  The  Cray  supercomputer  was  required  to 
perform  these  calculations.  Since  the  solution  of  the  original  pair  of  integral 
equations  had  required  most  of  the  computer  time,  this  procedure  for 
avoiding  numerical  differentiation  has  more  than  doubled  the  computing 
effort  on  this  problem. 

In  the  present  paper,  simple  resolutions  of  both  difficulties  are  pre¬ 
sented.  For  the  first,  the  peak  is  removed  by  a  procedure  suitable  for  arbitrary 
shapes,  and  the  variable  of  integration  is  changed  so  that  many  more  points 
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are  concentrated  in  the  neighborhood  of  the  peak,  without  sacrificing  either 
the  periodicity  or  the  interval-uniformity  condition  for  applying  the  MAQF. 
For  the  second,  the  variable  of  differentiation  is  changed  so  that  the  added 
masses  vary  nearly  linearly  with  that  parameter  when  the  gap  is  small.  In 
terms  of  that  variable,  the  errors  of  numerical  interpolation  and 
differentiation  are  much  less  than  with  the  original  one. 

The  aforementioned  change  of  differentiation  variable  is  suggested  by  a 
parametric  form  of  Herman's  result  [10]  for  the  successive  doublet  strengths 
in  the  circles,  which  was  also  used,  but  for  other  purposes,  by  Mitra  [11],  and 
for  a  special  case  which  will  be  described  subsequently,  by  Shail  [12].  By 
applying  this  form,  efficient  procedures  for  summing  the  doublet  series  for 
the  added  masses,  and  asymptotic  formulas  for  the  added  masses  and  their 
derivatives  at  small  gaps  were  developed.  These  are  new  results  in  the 
classical  theory  of  successive  images,  which  was  considered  to  have  been 
solved  completely  by  Hicks  [13,  14]  and  Herman  [10]  over  one  hundred  years 
ago.  Their  papers  on  the  motions  of  two  spheres  along  the  line  of  centers 
present  two  different  closed-form  solutions  for  the  positions  and  strengths  of 
the  successive  doublets.  Their  results  for  the  positions  are  also  valid  for  the 
motions  of  two  circular  cylinders,  and,  by  a  slight  modification  of  their 
analysis,  the  strengths  of  the  corresponding  doublets  can  also  be  obtained. 
However,  their  derivations  are  difficult  to  comprehend,  especially  that  of 
Hicks,  partly  because  they  did  not  take  full  advantage  of  the  theory  of 
continued  fractions  and  recurring  series  which  was  well  known  at  the  time  of 
their  papers.  Although  a  clearer  derivation,  based  on  these  algebraic  theories, 
had  been  devised  by  one  of  the  present  authors  several  years  ago,  it  was  not 
considered  worth  publishing  at  the  time;  but  now  that  we  have  significantly 
extended  the  classical  theory,  it  appears  appropriate  to  include  these  devel¬ 
opments  as  part  of  a  new  derivation  of  the  successive-image  formulas.  Guo 
and  Chwang  [2]  also  show  Herman's  solution  [8],  but  without  derivation  or 
attribution  or  any  new  results  in  the  successive-image  theory  of  circles  or 
spheres  moving  along  their  line  of  centers.  They  also  claim,  that  a  "closed- 
form"  solution  for  the  added  masses  "will  be  derived",  but  that  did  not  appear 
and,  probably,  does  not  exist. 
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The  plan  of  the  present  paper  is  as  follows.  In  Part  I,  a  new  derivation 
of  Herman's  formulas  for  the  added-masses  by  the  method  of  successive 
images  will  be  presented,  including  a  new,  efficient,  and  accurate  procedure 
for  computing  the  added  masses,  which  also  yields  their  asymptotic  formulas 
when  the  gap  between  the  bodies  is  small.  Since  the  added  masses  are 
obtained  from  an  infinite  series  of  doublet  strengths,  and  this  series  must  be 
differentiated  term-by-term,  it  is  important  to  investigate  its  uniform  con¬ 
vergence.  This  was  undertaken  by  Guo  and  Chwang  [2,3],  but  their  treatment 
is  incomplete,  proving  only  a  necessary  condition  for  ordinary  convergence. 
Hence  proofs  of  uniform  convergence  will  be  presented. 

In  Part  n,  the  method  of  integral  equations  for  obtaining  the  added 
masses  will  be  treated.  Procedures  for  reducing  the  peaks  of  the  kernel  and 
discretizing  the  integral  equations  by  means  of  the  MAQF,  but  with 
nonuniform  intervals  which  concentrate  points  in  the  neighborhood  of  the 
peaks,  will  be  presented.  Lastly,  a  procedure  for  obtaining  accurate  results  for 
the  added-mass  derivatives  by  numerical  differentiation  will  be  described. 
Added  masses,  their  derivatives,  and  interaction  forces  for  various  cases  will 
be  computed  and  compared  with  the  ’exact'  results  from  Part  I. 
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Part  I:  Solution  for  a  Pair  of  Circles  or  Spheres  by  Method  of  Successive 

Images 

1.  Recurrence  formulas  for  doublet  strengths:  added  masses 

Circles  of  radii  a  and  b  are  moving  with  velocities  Ui  and  U2  along 
their  line  of  centers  which  will  be  taken  at  the  x-axis.  Denote  the  centers  of 
the  circles  by  A  and  B,  and  their  x-coordinates  by  xi  and  X2.  Put  c  =  X2  -  xi;  see 
Fig.  1. 


The  fluid  is  assumed  to  be  invisdd  and  incompressible  and  the  flow  ir- 
rotational.  The  velocity  potential  for  the  fluid,  which  is  assumed  to  be  at  rest 
at  infinity,  and  contain  no  other  stationary  or  moving  boundaries,  may  be 
written  in  the  form 


O  =  +  U2<)>2 


(1) 


where  cj)-,  is  the  velocity  potential  when  Ui  =  1  and  U2  =  0,  and  <t>2  that  when 
U2  =  1  and  Ui  =  0.  In  terms  of  separate  polar  coordinate  systems  (ri,0i)  and 
(r2,02)  with  centers  at  points  A  and  B  and  0i  and  02  measured  from  the  x-axis, 
the  boundary  conditions  are 


9ri 


a  -  «*  81. 


dr2  /2=b 


— k  =  0; 


9ri  J3 


fd<b  2\ 

=  '  (ar2)>  = 


COS  02 


(2) 


We  shall  apply  the  well-known  method  of  successive  doublet  images  at 
alternate  inverse  points  in  the  circles  to  determine  added  masses  and  forces. 
For  this  purpose,  we  need  to  derive  analytical  expressions  only  for  the  case  Ui 
=  1,  U2  =  0,  since  that  for  Ui  =  0,  U2  =  1  could  be  obtained  by  permuting 
indices.  For  the  doublets  within  the  circle  about  A,  let  d2n,  52n,  n  =  / 

denote  the  distances  from  B  and  their  strengths,  respectively.  Also,  for  the 
doublets  within  the  circle  about  B,  let  d2n-l/  S2n-1/  n  =  1,  2,  3,...  denote  the  dis¬ 
tances  from  A  and  their  strengths.  Then  we  have,  successively, 

do  =  C,  60  =  a2,  &2n  =  -  §2n-l  52n'1  =  *  ^2n~2  ^5^2^ 


which  gives 
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X.  X  (  ab  N2  Op(ab)^  _ 

8211  ^  Ivd2n-2d2n-l  j  (dodid2...d2n-l)2  R  '  '  '' 

Then  we  obtain  from  (3) 

(ab)211 


(4) 


82n-l  =  - 


(dodid2...d2n-2)2 


(5) 


Here  the  d’s  are  successively  given  by  the  inversion  formulas 


d2n  =  c- 


d2n-l 


d2n-l  =  c 


b2 

d2n-2 


(6) 


The  doublet  strengths  can  then  be  computed  successively  from  (3). 


When  circle  B  is  moving  with  unit  velocity  in  the  positive  x  -  direction 

*  *  * 

and  circle  A  is  at  rest,  new  sets  of  d's  and  8's  are  generated.  Let  d2n_|,  8211-1 '  n  = 

1,2,3,...  denote  distances  from  A  and  the  doublet  strengths  at  those  locations 

»  * 

within  circle  B,  and  d2n,  8^,  n  =  0,1,2,...  the  distances  from  B  and  doublet 

strengths  within  circle  A.  Then  (3)  and  (6)  are  modified,  by  interchanging  a 
and  b,  to  become 

d^  =  c,  S^  =  b2,  4-1  -  -  4-2  <-T“>2  (7a) 


*2n-l 


A2n-2 


*  b2 

d2n  =  C--^ 
d2n-l 


*  a2 

^n-l  =  c'  * 

d2n-2 


(7b) 


Similarly, 


Sgtab)21' 


52n  = 


*  *  * 


(do  d,  dj-.d^ 


c>*  (ab)211 

Y2n-1  =  ’  *  *  *  * 

(do  dl  d2*"d2n-2^2 


(8) 
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Although  the  recurrence  formulas  are  identical  with  those  of  (3)  and  (6),  the 
change  in  the  initial  values  yields  different  sets  of  numbers  for  the  d's  and  8's. 

As  is  shown  in  [1],  three  added  masses.  An,  A22  and  A12,  occur  in  the 
equation  of  motion  of  a  pair  of  bodies  along  their  line  of  centroids  when  the 
fluid  is  at  rest  at  infinity.  Here  An  and  A22  can  be  obtained  from  the  general¬ 
ized  Taylor  formulas  of  Landweber  [4]  for  the  three  dimensional  case,  or  of 
Landweber  and  Yih  [5]  for  the  two-  and  three-dimensional  cases.  These  give 
expressions  for  the  added  masses  of  a  moving  body,  when  all  other  bodies  and 
boundaries  are  at  rest,  in  terms  of  the  sum  (or  integrals)  of  the  source  mo¬ 
ments  and  doublet  strengths  within  or  on  the  surface  of  the  moving  body. 
Reference  [6],  by  Landweber  and  Chwang,  gives  a  similar  expression  for  A12  in 
terms  of  the  sources  in  or  on  a  stationary  body  induced  by  the  moving  one. 
Guo  and  Chwang  [2]  noted  that  the  contribution  from  the  doublet  strengths 
was  not  shown  explicitly  in  [6]  and  included  it  in  their  formulation.  That  was 
an  oversight  in  [6],  although  the  source-moment  formula  could  be  considered 
to  include  the  doublet  strengths  since  the  strength  of  a  doublet  is  defined  by 
the  moment  of  its  source-sink  pair.  The  resulting  expressions  for  the  added 
masses  for  the  two  circular  cylinders  are  then 

An  =  Tcp  (2 Zq  52n  -  a2)  =  Ttp  (a2  +  2 820)  (9) 

A22  =  Ttp  (2Xq  8*n  -  b2)  =  Ttp  (b2  +  2Z~  8^)  (10) 

A12  =  2itp  L7  8*n  l  =  27tp  82n-l  (ID 

Here  p  is  the  two-dimensional  (2-D)  mass  density  of  the  fluid. 

If  the  bodies  are  spheres  instead  of  circular  cylinders,  then  the  positions 
of  the  successive  doublet  images  are  also  given  by  (6)  and  (7),  but  the  strengths 
of  the  doublets  are  altered.  Instead  of  (3),  (4),  (5),  (7)  and  (8),  we  would  have 

do  =  c  80  =  2  a3  82n  =  -  S2n-1  (^)3  52n-l  =  "  &2n-l  ^d^2^  ^ 
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Then 


<4') 


Mab)"*1  lr  (ab)n 


^=-5o[dS^^3=-2fd«I: 


)33 


(50 


and 


d;=c  i4*  4 = -  4.,  <^->3  4-,  -  -  4.2  <^->3  <7'> 


♦  * 


(ab)n 


^n-l  u2n-2 

*  1 r  (ab)n 


52n  =  So  [~— ^r~~»  33  Szn.!  =  -  i  t  7 . 33 

dQ  d^-.d^.i  dD  d1...d2n_1 


(80 


For  the  spheres,  the  three-dimensional  form  of  the  added-mass  for¬ 
mula  must  be  used.  In  terms  of  the  masses  of  displaced  fluid  of  the  spheres, 
4  4 

Mi  =  g  rcpa3  and  M2  =  ^  rcpb3,  where  p  is  the  mass  density  of  the  fluid,  we 
obtain 


An  =  4rcp  Zg  82n  -  Ml  =  4jcp(l7  &2n  +  \  a3)  =  (1  +  Jj  82n) 


(90 


M2  ,  6  oo  *  , 

A22  =_2~(1  +PZ1  O 

(100 

A12  =  4ltp  S2n-1  =  47tp  ^n-l 

(110 

Because  the  subsequent  derivations  of  the  d’s  and  5's  for  the  circular 
cylinders  and  the  spheres  are  very  similar,  only  those  for  the  former  will  be 
presented,  and  the  corresponding  results  for  spheres  will  be  collected  in  a  sep¬ 
arate  section. 
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2.  Forces;  added-mass  derivatives 


Once  the  added  masses  are  known,  the  interaction  forces  Fi  and  F2  on 
the  bodies  can  be  obtained  from  Lagrange's  form  of  the  equations  of  motion 


d  9T  9T  _ 

* 

d  9T  3T 
dt  9U2  9x2  ^ 

as  is  shown  in  Ref.  [1].  Here  T  denotes  the  kinetic  energy  of  the  fluid, 

2T  =  AnUi  +  2A12  U1U2  +  A22U2 


(12) 

(13) 


(14) 


Expressions  for  Fi  and  F2  given  in  [1]  reduce,  for  the  special  case  that  will  be 
considered  in  the  present  paper,  that  Ui  and  U2  are  constants,  to 

1  2  •  1  •  ^ 

Fi  =  2  AnU|  -  A11U1U2  -  (A12  +  2  A22)  U2  (15) 


12  1 
F2  =  (2  A11  +  A12)  Uj  +  A22  U1U2  -  2  A22  Uj  (16) 

where  the  dot  over  a  symbol  indicates  differentiation  with  respect  to  c.  In  the 
more  general  formulas  for  Fi  and  F2,  both  the  added  masses  and  their 
derivatives  appear. 

In  any  case,  the  added-mass  derivatives  play  an  important  role.  These, 
obtained  by  differentiating  (9),  (10)  and  (11),  are  given  by 

* 

An  =  27tpz7  52n  A22  =  27CpE7  8^ 


(17) 

•  00  •  00  •  * 

a12  =  2jtp£|  82n-l  =  2ftp£|  92n.i 
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and  recurrence  formulas  for  obtaining  successive  values  of  the  5’s  can  be 
derived  by  differentiating  (3),  (6),  (7)  and  (8).  This  gives  for  circles 


52n  =  (dfr-i  <^2n‘1  '  ^2n'1^  ^2n-l  =  ^  d2n-2  "  ^2n'2)'  =  0 


A2n-1 


•"20-2 


d2n  =  1  +  (^^')2  d2n-l/  d2n-l  =  1  +  (^~^)2  d2n-2,  do=l 


(18) 


Similarly  we  obtain  for  spheres 


j.  a3  r3§2n-l  ^ 

52n  =  _^1  d^T'520’1 


•*2n-l 


•  b3  ( 382n-2  •  ^ 

&2n-l  =  ~3  dfT ' 52n'2  5°  = 0  (18,) 

u2n-2 


d2n-2  V 


with  the  same  equation  as  in  (18)  for  the  d's.  The  formulas  for  the  8*  and  d* 
are  the  same,  except  that  a  and  b  are  interchanged. 


3.  Solutions  for  the  locations  and  strengths  of  the  doublets 
We  have,  by  (6), 

b2 


b2 


do-C,  di  -  C  -  c  d2  -  c  '  d|  “  C  '  d3-C-‘  a2 


C  - ' 


c  - 


c-bVc 


and  in  general,  in  the  usual  notation  for  a  continued  fraction. 


b2  a2  b2  b2 

2n‘1  ~  C  Cl-  C2-  C3-  C2n-1 


a2  b2  b2 
2n  ~  C  ’  Cl-  C2-  C2n 


(19) 


where  the  subscripts  on  the  c’s  serve  only  to  count  the  number  of  fractions, 
i.e.  ci  =  C2  =  C3  =  ...  =  c. 


Put  d2n-l  =  c  -  n  1  d2n  =  C  -  (20) 

q2n-l  q2n 

Since  the  p's  and  q's  with  odd  and  even  indices  are  defined  by  different  con¬ 
tinued  fractions  in  (19),  these  do  not  satisfy  the  usual  recurrence  relation  for 


10 


successive  convergents  of  continued-fraction  theory.  An  alternative,  simple 
recurrence  formula  is  satisfied,  however,  and  this  will  now  be  derived. 

We  have,  by  (19), 

P2n-1_  b2  =  b2(cq2n-3  -  P2n-3) 
q2n-l  a2  (c2-a2)q2n-3  -  cp2n-3 

C"  _E2n^ 
q2n-3 

Then  we  may  put 

P2n-1  =  b^cq^  -  P2n-3)/  q2n-l  =  (c2-a2)q2n-3  *  cp2n-3 

Then 

q2n-3  =  ^  (P2n-1  +  b2p2n-3) 


and  substituting  this  and  the  corresponding  expression  for  q2n-i  into  the  sec¬ 
ond  of  the  previous  equations,  we  obtain  the  recurrence  formula 


P2n+i  =  a2p2n-i 

& 

CO. 

1 

(21) 

where 

a2  =  c2  -  a2  -  b2. 

•8 

II 

a. 

(22) 

Similarly,  eliminating  p2n-i  and  p2n-3  from  the  same  pair  of  equations,  yields 
the  same  recurrence  formula  for  the  q's, 

q2n+l  =  Ct2q2n-1  -  P4q2n-3  (23) 

Since  the  coefficients  a2  and  (34  are  symmetric  in  a  and  b,  this  recurrence  for¬ 
mula  is  also  applicable  to  the  p's  and  q's  with  even  indices;  i.e.  we  have 


P2n  =  a2P2n-2  -  P4P2n-4,  q2n  =  a2q2n-2  -  P4q2n-4  (24) 

With  the  initial  values 

p0  =  0  q0  =  1  P2  =  ca2  q2  =  c2  -  b2 
pi  =  b2  qi  =  c  p3  =  b^c2-!?2)  q3  =  ca2  (25) 


11 


the  values  of  d2n  and  d2n+i  can  be  successively  calculated  from  the  recurrence 
formula  (24),  although  the  formulas  given  in  (6)  are  more  convenient  for  this 
purpose.  By  means  of  the  former,  however,  exact  closed-form  expressions  for 
d2n  and  d2n+i  can  be  derived. 

It  is  shown,  in  the  theory  of  recurring  series,  that  a  sequence  such  as 
{p2n+ll  which  is  linearly  related,  with  constant  coefficients,  to  the  two  previ¬ 
ous  members  of  the  sequence,  has  a  generating  function  given  by 


Pi  +  (P3  -  a2Pi)x 
1  -  a2x  +  P4x2 


oo 


I  P2n+1  Xn 
n=0 


This  is  readily  verified  by  multiplying  the  equation  by  l-a2x  +  p4x2.  Then,  by 
(25), 


b2(l  +  a2x) 

r - 3 IP2n+lXn 

1  -  a2x  +  J34x2  n=o 

Similarly  we  obtain  for  the  other  three  sequences 


l-a2x  +  P<x2'n50q2’",!‘n 


ca2x 


=  I  P2n  Xn 


l-a2x+p4x2  n=i 
l+a2x 

i  2  R4  2  —  ^  ^l2n 

1  -  a2x  +  P4x2  n=o 
The  quadratic  denominator  may  be  written  as 

1  -  a2x  +  p4x2  =  (1  -  ^ix)(l-^2x) 

where 


A.i,2  =  2  (®2  +®i)/  cs\=X\-Xi  =  (a4 - 4p4)1/2 


Here  Xif2  are  real  since 


(26) 

(27) 

(28) 

(29) 

(30) 

(31) 
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and  then 


c2  >  (a+b)2 
c2  -  a2  -  b2  =  a2  >  2P2. 


Hence  writing 


1  1  ,  X2  v  1  ~  ,n  n 

- = - ( - - - )  =  —  Z  an+ixn,  an  =  X.  -  Xo 

(1-Xix)(l-X2x)  Xi-X2  1-Xix  1-X2x  01  o  ^  ^ 


and 


_ x _ 

(1-Xix)(l-X2x) 


1  00 

=  — Z-  GnXn 
Oi  1 


we  obtain  from  (26,  27,  28,  29,  30) 


P2n+1  =  —  (On+l  +  a2On), 
01 


q2n+l  * 


CC?n+l 

ai 


ca2dn 

P2n  =  /  q2n  = 

Ol 


<3n+l  +  a2On 

o\ 


(32) 

(33) 


Then,  by  (20), 


d2n  = 


CtTn+1 

O’n+l  +  a^°n 


d2n+l  = 


(c2-b2)orH-i  -  P4an 

COh+l 


Q’n+2+a^<7n+l 

COn+1 


(34) 


since  c2  -  b2  =  a2  +  a2  and  by  (31), 

X:2  =  a2Xi-p4,  i  =  1,2 

This  gives  the  locations  of  the  successive  doublets. 

In  order  to  find  their  strengths,  we  see  from  (4)  and  (5)  that  the  prod¬ 
ucts  of  the  d's  are  required.  From  (34),  we  obtain 


d2n-2  d2n-i 


Pn+1  +  a2gn 

On  +  a2on-i 


Hence 
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dodid2...d2n-l  = 


On+1  +  a^n 
Ol 


Also,  by  (34), 


(dodi...d2n-3)d2n-2  = 


gn+a^n-l  #  CCq  _  CGn 
Ol  On+a^On-l  ai 


Hence  (4)  and  (5)  give 


&2n  = 


L<Jn+l+a2Oi 


_‘2 

/ 

n_ 


&2n-l  =  ‘ 


L  COn  J 


(35) 


The  corresponding  results  for  the  case  that  Ui  =  0  and  U2  =  1  can  be  ob- 

♦  *  * 
tained  from  (34)  and  (35)  by  substituting  and  d2n+1  for  d2n  and  d2n+i/  8^ 

* 

and  52n_|  for  82n  and  82n-l/  and  interchanging  a  and  b,  i.e. 


*  cqn+i 

2n  On+l+l^On 


^2n+l  - 


On+2  +b2CTn+l 
con+l 


n  =  0,1,2,... 


(36) 


b^oi 

2 

* 

£ 

.gn+l+t^Gn. 

°2n-l  ~  " 

.  OTn  . 

(37) 


Results  (34)  and  (36)  for  the  locations  of  the  successive  doublets  were  first 
derived  by  Herman  [10]  for  a  pair  of  spheres.  With  a  slight  modification  of  his 
derivation,  one  can  obtain  the  results  (35)  and  (37)  for  the  strengths  of  the 
doublets  for  a  pair  of  circular  cylinders.  Here  only  the  derivation  of  the 
Herman  formulas  is  new.  These  are  also  presented  by  Guo  and  Chwang  [2], 
but  without  attribution  or  derivation. 


4.  Solutions  for  added-mass  coefficients 

The  nondimensional  added-mass  coefficients  are  defined  by 

ki=An/M2,  k2  =  A22/M2,  ki2  =  A12/M2;  M2  =  rcpb2  (38) 
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Then,  by  (9),  (10),  (11),  (35)  and  (37),  the  coefficients  are  given  by 


3.2  2  oo 

kl  =  g2  (1  +  ^2  ^1  ^2n)  =  ^2 


oo  ,  B^Oi  2 

1+2 2,  (— i£=— V~) 

L  an+i  +  a2an  . 


k2  =  l+^E752n  =  l+2i:7 


(  >2 


an+i+bZanJ 


,  2  _oo  ™  2  OO 

kl2  =  ^2i1  02n-l  =  -pii 


(^Hqi 
L  CCJn  J 


In  order  to  avoid  very  large  numbers  in  computing  ki,  k2  and  ki2/  we 
the  nondimensional  quantities 

Ti  =  a2/(32  a;  =  an/P 20  a^i/P2  ^  =  A.2/P2 

Then  (39),  (40),  and  (41)  become,  in  terms  of  xn  = 

2 

ki  =  £2  [1  +  2b2Z~(bxn+l  +  axn)'2] 


k2  =  1  +  2a2E^(axn+i  +  bxn)"2 


,  2a2 

k12  =  -'3"2:i  T 


-2 

n 


and  from  (31), 


a{  =  (ti2-4)1/2  X[  2  =  \  [ti  +(t|^-4)1/2] 

The  forms  in  (42)  suggest  the  substitution 

r|  =  2  cosh  £ 


which  yields 


=  2  sinh  £  =  e±<* 


a'  =  2  sinh  n£ 


(39) 

(40) 

(41) 

define 

(39a) 

(40a) 

(41a) 

(42) 

(43) 

(44) 
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Since,  hereafter,  only  the  nondimensional  forms  of  X],  X2  and  an  will  be  used, 
the  primes  on  these  quantities  will  be  omitted.  Then,  since 


sinh  n£ 
Xn  sinh  £ 


(45) 


by  (44),  the  added-mass  coefficients  may  be  written  in  the  third  form 


ki=p{l+2b2z7 


sinh2£ 


[b  sinh(n+l)£  +  a  sinh  n£]2 


} 


(39b) 


k2  =  1  +  2a2E7 


sinh2£ 


[a  sinh(n+l)£  +  b  sinh  n^]-1 


k12  =  -fr2®  i2(C)^T~g 


(40b) 

(41b) 


It  will  be  shown  that  the  parametric  forms,  (39b,  40b,  41b),  which  appear 
to  be  new  except  for  Shah's  result  described  below,  play  an  important  role  in 
summing  the  three  infinite  series  and  investigating  their  uniform 
convergence,  and  enable  the  asymptotic  behavior  of  the  added  masses  and 
their  derivatives  to  be  determined  when  the  gap  between  the  circles  is  very 
small. 


The  same  parameter  £  was  introduced  by  Mitra  [11]  and  Shail  [12]  in 
their  treatments  of  two-sphere  problems  by  the  method  of  successive  approx¬ 
imations  using  correction  potentials.  Mitra  solved  only  the  electrostatic 
condenser  problem.  Shail  applied  Mitra's  analytical  method  to  the 
hydrodynamic  problem  of  two  spheres  moving  along  their  line  of  centers,  but 
gave  specific  results  only  for  the  case  of  a  sphere  approaching  or  moving  away 
from  a  wall  along  its  normal.  His  result  for  k2  agrees  with  that  in  (75a)  when 
the  wall  is  treated  as  a  sphere  of  infinite  radius. 

We  shall  now  show  that,  even  when  £  is  small,  n  need  not  be  taken 
very  large  in  evaluating  the  infinite  series  (39b,  40b,  41b).  To  illustrate  the 
procedure  to  be  used,  let  us  apply  it  to  evaluate  numerically  the  series 
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1/n2.  We  shall  see  that  this  series  is  related  to  the  value  of  ki2  at  impact. 
It  is  well  known  that 


K2 

-g-  =  1.644934... 


(46) 


This  may  be  written  as 


L^  +  R*(N),  R*(N)  =  Zn+122 


Interpreting  the  remainder  R*(N)  as  the  trapezoidal  rule  for  the  integral 


R(N)  = 


oo 


dn  _ 
n2  - 


we  obtain  the  values  for  the  error  in  truncating  the  series  at  N  shown  in 
Table  1. 


Thus  the  sums  from  1  to  N  have  errors  of  one  in  the  first,  second,  and 
third  decimals  for  N  =  10,  100  and  1000  respectively,  but  essentially  agree  to 
four  decimals  at  N  =  10,  and  to  at  least  six  decimals  at  N  =  100  and  1000  when 
corrected  by  the  "remainder"  given  by  the  trapezoidal  integral.  Since  the 
added-mass  series  converge  faster  than  the  series  (46),  as  will  be  shown,  this 
suggests  that  adequate  accuracy  could  be  obtained  with  N  =  0(10),  provided 
the  n-th  terms  of  the  three  added-mass  series  could  be  integrated  into  a  closed 
form.  That  this  is  the  case  will  now  be  shown. 


Let  us  first  consider  ki2,  which,  by  (41b)  may  be  written  as 

2a2  n  sinh2C  2a2  °°  sinh2£  / 

«+r’2(n)'  Ri2(n,=  cT  <lS^dn  (47) 

N+2 


But 


2a2  °°  2a2  sinh2C  i 

Rl2(N)  =  -  —  sinh2C  coth  n£  L  i  =  — - ~  exp  [-(N  +  (48) 

c2C  ^+2  c2C  sinh(N+k 
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As  £  approaches  zero,  we  obtain  from  (48)  the  limiting  value 

,  .  2a2  oo  1  JC2a2  a2  / 

ki2(0)  =  -  ^2  =  "  ~^T  =  -  3.289868-  p  s  =  a+b  (49) 

The  limiting  value  of  [9ki2(C)/^C]c=0  can  also  be  obtained  from  (48).  Since 
each  term  of  the  truncated  series  is  an  even  function  of  £,  its  derivative  is 
zero  at  £  =  0.  Also,  as  £  — »  0  for  a  fixed  N, 

2a2  _  l  2a2  1 

—  sinh2C  [coth(N  +  2)  C  -  U  -  ( - I '  0 

c2C  ^  N  +  ^ 

Hence,  indicating  differentiation  with  respect  to  C,  hereafter  by  a  prime,  we 
have 


k12(0)  = 


2a2 

s2 


and 


2a2  7c2 

ki2(0--p-(-g-0 

Next  consider  k2(Q,  written  in  the  form 


k2(0  =  1  +  2, 


N 


2a2  sinh2^ 


[a  sinh(n+l)£  +  b  sinh  n£]2 


+  R2(N,C) 


(50) 


(51) 


(52) 


where 


00 

R2  =  j 


_ 2a2sinh2£dn _ 

[a  sinh(n+l)£+b  sinh  n£]2 


2a  sinh  £  sinh  n£ 

£[a  sinh(n+l)£  +  b  sinh  n£] 


At  the  upper  limit,  we  have 


sinh(n+l)£  e(n+l)C 
sinh  n£  en£ 


18 


Then 


2a  sinh  £  2a  sinh  £  sinh(N+^)£ 

2  "  C(aeC+  b)  £[a  sinh(N+  |)C  +  b  sinh(N+  \)Q 

2a2sinh2C’exp[-(N+ 

= — ; - r -  (53) 

£(aeC+b)[a(cosh  C,  +  coth(N+2)C  sinhQ  +  b] 


Here  also  only  the  term  from  the  upper  limit  of  the  integral  contributes  to  the 
coefficient  of  £  in  the  Taylor  expansion  of  k2  about  £  =  0,  since  the  other  terms 
are  even  functions.  As  £  — »  0,  we  obtain 


R 2(0  - 


2a 

s+a£ 


2a  a 

-t*1 -?o 


(54) 


Then 


We  then  have 


k2(0  -  k2(0)  -  K  C 
& 


(55) 


(56) 


where,  by  (52)  with  N  =  °o  and  £  — >  0,  we  obtain  the  convergent  series 


k2(0)  =  1  +  17 


2a2 

(sn+a)2 


(57) 


Similarly,  we  find 


2b2sinh2^ 


2a2 


k'®  '  *  11  +  [b  sinh(n+l)£  +  a  sinh  n£j2^  "  k>®  "ST  C 


(58) 


where 


n  T-  2b2  , 
[1  +  Z*  (sn+b)2 


(59) 
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The  linear  variation  of  k-j  and  k12  with  £  is  shown  graphically  in  Fig.  2  for  the 
case  a  =  b,  for  which  (58)  and  (59)  become 

ki(0)  =  k2(0)  =  1+  2^7  (2n^1)2  =  X  ’  1  =  1-467401-  (60) 

The  results  for  ki(0),  k2(0)  and  ki2(0),  corresponding  to  c  =  a+b,  were  given  by 
Guo  and  Chwang  [2].  The  expressions  for  k{(0),  k£(0)  and  k12(0),  however,  are 

believed  to  be  new.  We  see  that  the  three  added-mass  coefficients  have  slopes 
of  the  same  magnitude  at  £=0.  The  remainder  Ri(Q,  obtained  from  (53)  by 
interchanging  a  and  b  and  multiplying  by  aVb2,  becomes 

2a2sinh2C»exp[-(N+ibC] 

Ri(Q  =  — - - r -  (6i) 

C(be£+a)[b(cosh£+coth(N+2)C  sinhQ+a] 


5.  Asymptotic  formulas  for  small  gaps 

The  principal  mathematical  parameter  in  the  foregoing  equations  is  £, 
which  was  defined  in  (43).  The  principal  physical  parameter  is  the  gap  g  be¬ 
tween  the  circles,  given  by  c  =  s+g.  When  g/b  and  a/b  are  given,  then  r\  and  Oi 
=  Cq2-4)l /  2  are  also  known,  and  C,  can  be  found  by  (43)  or  (44)  by  inverting 
either  the  hyperbolic  cosine  or  sine.  When  g  is  relatively  small,  however,  the 
relation  between  these  parameters  becomes  very  simple  and  useful.  This 
relation  will  now  be  derived. 

We  define  the  nondimensional  quantities  y  and  (i  by 

yffi"  ^=7  mV-?  «2) 


Then,  we  obtain 


a.2  = 


and,  by  (43), 
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(63) 


cosh  C  =  1+  ^Y2  +  ^ H2'/* 


Since  cosh  C,  =  1  +  2sinh2  (C/2),  then 

C  =  2sinh‘1Z  =  2[z  -  \  Z3  +  J^Z5  Z<1  (64) 

where,  by  (63), 

Z  =  2  (1  +  \  M-2)^) 1/2  =  |  (1  +  |  vW-  nV  +  *“)/  M-V  =  if 


Substituting  the  series  for  Z  into  (64),  we  now  obtain  the  asymptotic  formula 
C  - y[i  - ^ d  -  3m-2)  72  +  6io(3-  10p2  - 5p4)  y1]  (65) 


Of  the  two  convergence  conditions,  Z  <  1  and  2g/s  <  1,  the  former  is 
more  stringent.  The  latter  restricts  g  to  the  mean  of  the  radii,  g  =  (a  +  b)/2. 
The  former  can  be  expressed  as 


g 

s 


< 


1  + 


4abN 


1/2 


-1 


in  which  the  right  member  has  a  maximum  value  of  y[l  -  1  when  a  =  b,  and  is 
zero  when  the  ratio  of  the  radii,  a/b,  is  zero  or  infinity.  Even  at  the 
maximum  value,  we  have  g/a  <  -  1)  =  0.828...  versus  1.0  from  the  other 

condition.  Yet,  as  we  shall  see,  (65)  yields  a  useful  approximation  at  g/a  =  1. 
Although  the  infinite  series,  of  which  (65)  is  a  truncation,  diverges,  that  series 
belongs  to  the  class  of  divergent  series  in  which  the  error  in  using  a  partial 
sum  is  less  than  the  next  term  of  the  series.  Thus  the  partial  sum  is  useful  if 
the  next  term  is  very  small. 


According  to  (62)  and  (65),  C  is  proportional  to  the  square  root  of  the  gap 
g  when  the  gap  is  very  small.  For  example,  (58)  gives 

ki(C)  =  k!(0)-^(^)1/2  (66) 
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Then 


dlq  a2  2s  1/2 
dg  *  s2  (p2g) 


— >  oo  as  g  0 


in  agreement  with  the  known  property  that  the  slopes  of  the  added-mass 
coefficients  versus  g  become  infinite  as  the  gap  approaches  zero.  That  the 
asymptotic  slope  varies  inversely  as  the  square  root  of  the  gap,  however,  is 
believed  to  be  a  new  result. 


For  small  values  of  the  gap,  the  asymptotic  formula  (65)  yields  highly 
accurate  values  of  C-  For  example,  with  a  =  b,  we  have  p.  =  1  /2  and  y  =  2 V  g/a. 
Then,  at  g/a  =  0.25,  (y  =  1),  (65)  gives  C  =  0.989876  and  sinh  C/2  =  0.515394 
versus  the  exact  value  from  (64),  0.515388.  The  agreement  is  even  better  at 
smaller  gaps.  At  g/a  =  1,  (65)  gives  C  =  1.92604  and  sinh  C/2  =  1.1189  versus  the 
exact  value  1.1180.  Thus  the  asymptotic  formula  for  C(y)  is  sufficiently 
accurate  for  values  of  g/a  <  1. 

Equations  (62)  through  (65)  are  also  valid  for  two  spheres.  The 
asymptotic  law  for  that  case,  given  in  (78),  differs  from  that  for  the  circles  in 
(66). 

6.  Derivatives  of  the  added-mass  coefficients 

Recurrence  formulas  for  calculating  the  derivatives  of  the  added 
masses  were  given  in  (18).  Since  the  accuracy  may  be  improved  by  applying 
(39a),  (40a)  and  (41a),  expressions  for  the  derivatives  will  now  be  derived  on 
the  basis  of  these  equations.  Derivatives  with  respect  to  C  will  now  be  indi¬ 
cated  by  the  prime  of  a  quantity;  e.g.,  dx/ dC  =  x'.  Then,  from  (39a,  40a  and  41a), 
we  obtain 


k|  =  -  4a2  (bxn+i  +  axn)’3  (bx^  +  ax£)  (67) 

k£  =  -  4a2  (axn+i  +  bxn)'3  (ax^  +  bt^  (68) 

cUt) 

k/2  =  4a2  (ctn)'3  (c'tn+l  +  ex'),  c'  -  —  sinh  C  (69) 
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in  which  c'  is  obtained  by  differentiating  a2  =  2|32  cosh  Also,  by  (45),  we 
have 


=  csch  £  (n  cosh  n£  -  rn  cosh  Q  (70) 

with  which  the  k'’s  can  be  computed. 

The  trapezoidal-integral  procedure  for  evaluating  the  remainder  of  a 
truncated  infinite  series  may  also  be  applied  to  the  series  of  derivatives  of  the 
added-mass  coefficients.  To  apply  the  method,  it  was  necessary  that  the  n-th 
term  of  the  series  be  integrable  in  closed  form  with  respect  to  n,  as  was  shown 
to  be  the  case  for  the  parametric  forms  of  the  coefficients.  Hence  the  n-th 
term  of  the  series  of  their  derivatives  must  also  be  so  integrable  since  the 
order  of  performing  differentiation  with  respect  to  £  and  integration  with 
respect  to  n  may  be  reversed;  i.e.,  the  integral  of  the  derivative  is  given  by  the 
derivative  of  the  integral  that  has  already  been  obtained  for  the  parametric 
form  of  the  added-mass  coefficients. 

The  method  will  now  be  demonstrated  by  evaluating  1,^(0 ■  We  obtain 

by  (49) 


,  .  n  d  (  sinh2£  ^  d  sinh2£  coth  n£ 
1  d£  (sinh2n£ J  d£  C 


N+; 


.  „n2,  ,  „  ,  sinh2^ 

=  ~2  (coth£  -  n  coth  nQ  -  — — 

T_  S 


•  {(2  C  coth  C  -  1)  [1  -  coth  (N  +  |)C1  +  (N  +  |)  C  csch2(N  +  |)c} 

7.  Convergence  of  the  series  for  the  added  masses  and  their  derivatives 

In  summing  the  series  for  the  added-mass  coefficients,  it  was  observed 
that  the  terms  of  the  series  were  even  functions  of  £ ,  so  that  the  derivative  of 
each  term  with  respect  to  £  is  zero  at  £  =  0.  Thus,  if  the  series  of  derivatives 
were  uniformly  convergent,  its  sum  would  also  be  zero  at  £  =  0,  contrary  to 
the  results  given  in  (51),  (55)  and  (58).  It  is  important,  then  to  investigate  the 
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convergence  of  the  series  for  the  added-mass  coefficients  in  (39a),  (40a)  and  for 
their  derivatives  in  (67),  (68)  and  (69). 

This  topic  was  treated  by  Guo  and  Chwang  [2],  but  their  treatment  was 
flawed  by  several  errors.  One  is  that  they  refer  to  proving  the  uniform  con¬ 
vergence  of  the  sequence  of  the  terms  of  a  series  {an},  instead  of  the  sequence 
of  partial  sums  {X1}  ar}.  A  second  is  that  the  simple  ratio  test  an+i/an  <  1  was 

used  to  "prove"  uniform  convergence,  although  that  criterion  is  only  a 
necessary  condition  for  convergence;  for  example,  it  is  satisfied  by  the  diver- 

oo 

gent  series  X1  (1/n).  Thirdly,  no  attempt  was  made  to  prove  uniform 

convergence  of  the  original  series  in  the  closed  region  0  <  g  <  Its  conver¬ 
gence  at  g  =  0,  and  uniform  convergence  in  the  open  region  0  <  g  <  °°,  does 
not  assure  uniform  convergence  in  0  <.  g  <  ■*>  since  the  series  may  be  discon¬ 
tinuous  at  g  =  0.  Lastly,  although  they  are  aware  of  the  need  to  prove  uniform 
convergence  of  the  series  of  derivatives,  the  reason  for  undertaking  to  prove 
that  the  original  series  is  uniformly  convergent  is  not  mentioned.  Hence  a 
reconsideration  of  this  subject  appears  to  be  necessary. 

First  consider  the  series  for  ki2  in  (41b).  Since  sinh  n£  >  n  sinh  0  <  C, 

-2 

<  oo,  (as  is  seen  from  their  Taylor  expansions),  then  by  (45),  xn  >  n  and  xn  ^  1/n2 
for  0  <  £  <  oo.  Since  also  c  >  s,  then 

X~  (ccn)-2  <  X^  (sxn)-2  <  S'2  X^  1  /n2 


Application  of  the  Weierstrass  comparison  test  with  the  series  X~  l/n2  then 
shows  that  the  series  for  ki2  is  uniformly  convergent  in  0  <  £  < 

Next  consider  the  series  for  k2  in  (40a).  We  have 
(axn+i  +  bxn)'2  <  (bxn)'2  <  (bn)-2 

Hence,  by  the  same  comparison  test,  the  series  for  k2  is  uniformly  convergent 
in  0  <  £  <  o°.  Similarly,  by  interchanging  a  and  b,  we  can  show  that  the  series 
for  kj  has  the  same  property. 
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Since  the  terms  of  the  three  series  are  continuous  functions,  and  the 
series  are  uniformly  convergent  in  0  <  £  <  we  conclude  that  the  added- 
mass  coefficients,  as  given  by  their  series,  are  continuous  functions.  This  was 
not  obvious  since  the  terms  of  the  series  involve  sinh  n£  which  is  zero  when 
£  =  0,  no  matter  how  large  n  may  be,  but  becomes  infinite  as  n  approaches  in¬ 
finity  no  matter  how  small  £  may  be. 

To  investigate  the  convergence  of  the  series  for  k£,  we  consider  its  form 
given  in  (40b).  For  n  very  large,  we  have 

_ sinh  £ _  2e~n£sinh  C, 

a  sinh  (n+l)£  +  b  sinh  n£  aeC  +  b 


Then 


d  e*nC  sinh  L  .  2ne_2nC  sinh2C  ,  x 

(  aeC  +  b  )2  '  -  (aeUb)?  <71) 

tfol l  .41AUA  '*!' 

Put  p  =  e'2C,  and  let  N  be  a  large  value  of  n  at  which  (71)  is  valid.  Then  (71)  is 
seen  to  be  proportional  to  the  derivative  of  the  power  series  pn  which, 

together  with  all  its  derivatives,  is  uniformly  convergent  within  its  ’circle  of 
convergence’,  i.e.,  0<^c°°,  orpcl.  A  similar  proof  verifies  that  the  series 
for  k{  in  (39b)  is  uniformly  convergent  in  the  same  open  range. 


For  kj2  in  (41a),  the  n-th  term  of  the  series  is 

2CTl°n  2c'al  9rr  rr' 
d  oi  1  n  1  2°l°i 


+ 


COn  c2C„  C^o" 


n 


n 


c2a„ 


where  c'  is  given  in  (69).  Since  an  =  Xn^i  and  the  series  (l/xn)  has  already 

been  shown  to  converge  uniformly  in  the  closed  region,  we  see  that  the  series 
of  the  second  and  third  terms  also  converge  uniformly  in  that  region.  The 
first  term,  however,  becomes  asymptotically 
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Ics-^n  cosh  n£  4Cj 
c2  sinh2  n£  c2  ne 

proportional  to  that  in  (71),  for  which  uniform  convergence  in  the  open  re¬ 
gion  has  already  been  proved.  Hence  the  series  for  k12  is  also  uniformly 
convergent  in  the  open  region  0  <  £  < 

If  the  derivative  series  were  uniformly  convergent  in  the  closed  region, 
then  k{  (0)  (say)  would  have  been  zero  since  all  the  terms  of  the  series  are 

even  functions;  but  this  would  contradict  the  results  in  (58).  Hence  there  is  a 
discontinuity  in  slope  as  £  passes  through  zero  from  positive  to  (physically 
meaningless)  negative  values,  with  ki(£)  =  ki(-Q  and  a  sharp  peak  at  £  =  0. 


8.  Results  for  spheres 


Basic  iteration  formulas  for  the  locations  and  strengths  of  successive 

doublets  for  a  pair  of  spheres  were  given  in  eqs.  (3')  to  (11').  These  show  that 

the  only  changes  in  the  expressions  for  the  doublet  strengths  from  those  for  a 

pair  of  circles  are  a  factor  of  1/2  and  changes  from  squares  to  cubes,  correlated 

* 

with  the  changes  from  80  =  a2  to  a3/2  and  8Q  =  b2  to  b3/2.  The  formulas  for  the 
positions  of  the  doublets  remain  the  same.  Thus  we  obtain 


$2n  -  o  t: 


aboi 


-]Kk 


absinh£ 


2  bon+i+aon  ^  sinh(n+l)£+a  sinh  n£ 

1  abqi3  1  r  ab  sinh£  ,3  * 

02n-l  —  "  2  l  J  —  “  2  L  .  ,  r  *  ~  ^2n-l 
z  CGn  z  c  sinh  n(. 


(72) 

(73) 


The  expressions  for  S2n  and  82n.|Can  be  obtained  from  (72)  and  (73)  by 
interchanging  a  and  b. 

The  added-mass  coefficients  can  now  be  obtained  by  substitution  into 

4 

[9"],  [10']  and  [11 ']  and  dividing  by  M2  =  g  Ttpb3.  This  gives 
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(74) 


,  ,  _o°  _  _ b  sinh  ^ _ 3 

1  2b3  ^  ^  b  sinh(n+l)£  +  a  sinh  n£ 


1  ,  *  r  a  sinh  C  ,3, 

k2  =  2{l+3L1  [  ...  ^  7-  .  ~  ~]  } 

z  a  sinh(n+l)Q  +  b  sinh  nc. 


(75) 


,  3a3  _  -  sinh3  C 

kl2  -  -2c3  m>  ^3  -  sinh3  n? 


(76) 


The  general  terms  of  the  series  (74,  75,  76)  can  also  be  integrated  into 
closed-form  expressions,  so  that  the  technique  of  applying  the  integral  to 
evaluate  the  remainder  after  the  N-th  term  may  be  used  to  obtain  high 
accuracy  with  a  moderate  value  of  N.  For  (76),  the  integral  to  be  evaluated  is 
of  the  form 


J  sinhrxdx 

for  which  a  reduction  formula  for  positive  or  negative  integral  values  of  r  is 
given  in  tables  of  integrals.  In  the  present  case,  r  =  -  3,  and  the  reduction  for¬ 
mula  reduces  it  to  the  known  integral  for  r  =  -  1.  The  denominators  of  the 
terms  of  (74)  and  (75)  can  be  expressed  in  the  same  form  by  writing  that  in 
(75),  for  example,  in  the  form 

a  sinh  (n+1)  £  +  b  sinh  n£  =  (A  sinh  n£  +  B  cosh  nQ  =  VA2-B2  sinh  (n£+p) 

A  =  a  cosh  £  +  b,  B  =  a  sinh  £,  coth  p  =  A/B. 

The  actual  integrals  will  be  presented  here  only  for  (76)  since  we  are  mainly 
interested  in  the  numerical  results  for  pairs  of  cylinders  in  the  present  paper. 

As  in  (49),  we  write 


3a3  n  sinh3^ 

’  12  "  2c3  [Ll  sinh3n£ 


+  Ri2(N)], 


oo 


Rl2(N)  = 


sinh3^ 

sinh3n£ 


and  we  find 
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Ri2(N)  +  —  ^  -  [csch  (N  +  C,  coth  (N  +  ^)  £  +  In  tanh  £] 


(77) 


As  in  the  derivation  of  (51),  the  limiting  value  of  (dki2/dQ^=o  can  be  obtained 
from  (76)  and  (77).  Again  we  see  that  all  the  terms  of  the  truncated  series  and 
the  first  term  of  Ri2(N)  are  even  functions  of  £,  so  that  only  the  ln-term  can 
contribute  to  the  derivative  at  £=0.  Then,  for  small  values  of  £,  as  C,  ap¬ 
proaches  zero  for  a  fixed  N,  we  obtain 


3a3 

kl2©  -  k,2<0)  -  53  C2  In? 

(78) 

3a3  3a3 

12-2^  C  In  C- 233  Yin  T 

(78a) 

since  according  to  (65).  This  shows  that  k12  approaches  zero  as  C, 
approaches  zero.  The  derivative  with  respect  to  c,  however,  is  infinite  since, 
by  (69), 

asC-»0.  (79) 

Similarly  it  can  be  shown  that  the  derivatives  of  ki  and  k2  with  respect  to  £ 
are  zero,  and  with  respect  to  c  are  infinite.  Again,  this  implies  that  the 
spheres  would  never  meet. 

The  series  for  the  added-mass  coefficients,  given  in  (74),  (75)  and  (76), 
converge  for  all  values  of  £,  0  <  £  <  At  £  =  0,  we  obtain 

kl(0)=^{1+3Z7  [naX1)b]3)  (BO) 


k2(0)  "  2  {1  +  3E1  [(n  +l)a+nb]3} 


(81) 


ki2(0)  = 


3a3 

‘2s3 


1.2020569... 


3a3 

2s3 


(82) 
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When  a  =  b,  we  get 


ki(0)  =  k2(0)  =  \  [1  +  3E“  (2n^1)3]  =  0.5776997...  (83) 

The  results  in  (80)  to  (83)  were  given  by  Hicks  [13]. 

The  proofs  of  uniform  convergence  of  the  series  for  ki(£),  k2(Q  and 
ki2(C)  in  the  closed  region  0  <.  C  <  00  for  the  two  spheres  are  similar  to  those 
for  the  two  circles,  with  the  convergent  series  in  (82)  serving  for  the  compari¬ 
son  test.  The  proofs  of  the  uniform  convergence  of  the  series  for  dki/d£, 
dk2/d£  and  dkj2/ d£  in  the  open  region  0  <  £  <  °°  are  also  similar,  with  p  de¬ 
fined  as  p  =  e~3C,  the  proof  then  depending  on  the  uniform  convergence  of  a 
power  series  in  p,  and  its  derivatives  within  its  'circle  of  convergence.' 


Fart  II;  Solution  by  Means  of  Integral  Equations 
1.  Derivation  of  integral  equations 

The  developments  described  in  Part  1  were  undertaken  to  evaluate  the 
accuracy  of  the  integral-equation  procedures  described  in  References  [1],  [2],  [3], 
and  [6],  and  used  to  obtain  added  masses  for  body  pairs  of  various  shapes. 
Among  these  procedures,  the  two  that  most  affected  numerical  accuracy  are 
the  treatment  of  the  sharp  peaks  of  the  integrands  of  the  integral  equations 
when  the  gap  between  the  bodies  is  very  small,  and,  secondly,  the  procedure 
for  computing  the  derivatives  of  the  added  masses. 

The  unknown  functions  in  these  integral  equations  are  taken  to  be  the 
source  distributions  ai  and  a2  on  the  surfaces  of  bodies,  moving  with  veloci¬ 
ties  Ui  and  U2  along  the  line  joining  their  centroids.  When  Ui  =  1  and  U2  = 
0,  the  solutions  Oj  and  a2  of  the  pair  of  integral  equations,  derived  from  the 

associated  boundary  conditions,  give  the  added  masses  by  applying  the  gener¬ 
alized  Taylor  formulas  [4],  [5]  or  [6].  The  x-axis  is  taken  along  the  line  joining 
the  centroids  and  the  bodies  are  assumed  to  have  an  x-z  plane  of  symmetry  so 
that  they  can  move  without  rotation  along  the  x-axis;  see  Fig.  3. 
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For  the  two-dimensional  case,  the  velocity  potential  <|>(P)  at  a  point  P  ex¬ 
terior  to  both  bodies  or  on  their  bounding  contours,  is  given  by 

4>(P)  =  oi(Qi)  In  rpQi  dsQi  +  ^  02(Q2)  In  rpQ2  dsQ2  (84) 

in  which  Qi  and  Q2  denote  points  on  the  closed  contours  ci  and  C2  of  bodies  1 
and  2  respectively,  rpQi  and  rpQ2  are  distances  to  P  from  Qi  and  Q2,  and  dsQi 
and  dsQ2  are  elements  of  arc  along  ci  and  C2-  The  boundary  conditions  to  be 
satisfied  are 


80  ^  3x 

8npi  3npi 


3d)  8x 

— — =  U2 - 

3np2  3np2 


(85) 


where  npi  and  np2  denote  distances  along  the  outward  normal  to  ci  and  C2, 
and  the  point  P  lies  on  the  boundary  of  body  1  or  2.  These  give  the  pair  of  in¬ 
tegral  equations 


-  a 

tcoi(Pi)  +  9^  oi(Qi)  ^77ln  rpiQi  dsQi  + 


8npi 


*<2  °2(Q2)  ^ ln  rP]Q2  dSQ2  =  Ul 


(86) 


t  a 

7co2(P2)  +  9C2  02(Q2)  ln  rp2Q2  dsQ2  + 


3np2 


oi(Q,,dklnrP2Q,dSQ1=U23^ 


(87) 


Here  the  additional  terms  7tGi(Pi)  and  7ta2(P2)  appear  in  differentiating  the  in¬ 
tegrals  in  (84)  which  are  singular  when  Qi  passes  through  Pi,  or  Q2  passes 
through  P2.  In  (86)  and  (87),  however,  the  kernels 


K(Pi,Qi)  -  On  rp-jQ^,  K(P2,Q2)  -  (In  rp2Q2) 


Pi 


8n 


(88) 


P2 


become  indeterminate  and,  at  smooth  points  of  the  contours,  have  the  limit¬ 
ing  values 
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(89) 


Ki(Pi,Pi)  =  \  Ci(Pi),  K2(P2,P2)  =  \  C2(P2) 


where  Ci  and  C2  are  the  curvatures  of  contours  ci  and  c2;  see  ref.  [15]. 

When  ci  and/or  c2  are  not  defined  analytically,  there  would  be  a  loss  of 
accuracy  in  determining  Ci  and  C2.  This  can  be  avoided  by  taking  advantage 
of  the  well  known  property  of  the  integral  of  the  transpose  of  the  kernel, 

Ki(Qi,Pi)  dsi  =  7i,  i  =  1,2  (90) 


Then  (86)  and  (87)  may  be  written  as 


3  ,  a 

27tai(Pi)  +  [tfi(Qi)  ln  rpiQ!  -  o1  (pi>— rPlQl^  dsQl 


an 


PI 


an 


Q1 


f  a  ,  ax 

+  ^C2  °2(Q2)  rPlQ2  dsQ2  -  U1 


3n 


PI 


an 


(91) 


PI 


3  d 

27ta2(P2)  +  i>C2  [ct2(Q2)  —  In  rp2Q2  -  o2(P2)  ~  rp2Q2]  dsQ2 


3n 


P2 


3n, 


Q2 


a  dx 

+  ^  rP2Ql  dsQl  =  U2 


3n 


P2 


an 


(92) 


P2 


We  see  from  (91)  and  (92)  that  the  integrands  with  the  indeterminacies  van¬ 
ish  when  P  and  Q  coincide. 


The  integral  equations  for  the  three-dimensional  case  are  derived  in  a 
similar  manner.  The  potential  at  a  point  P  exterior  to  the  bodies  or  on  their 
bounding  surfaces  is 

4>(P)  =  -  Js!  <*i(Qi)  dSQi  -  a2(Q2)  dSQ2  (93) 

where  Si  and  S2  denote  the  surfaces  of  bodies  1  and  2,  and  dSQi  and  dSQ2 
their  area  elements.  The  boundary  conditions  are  the  same  as  in  (85),  and  ap¬ 
plying  these,  we  obtain  the  pair  of  integral  equations 
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2TO1(P,)  -  sJ  01(0,)  dSQ,  -  sJ  02(02)  dSQ2  =  U,  ^  (94) 

2JtO2(P2)-sJ02(Q2)^^dSQ2-s|O1(Q,)^^dSQ,=U2^  (95) 
In  this  case,  the  kernels 


K.aw»-4s5?  *-«  <96) 

remain  singular,  and  it  is  more  important  than  in  the  two-dimensional  case 
to  remove  the  singularity.  This  can  be  accomplished  by  applying  the  property 
of  the  transpose  of  the  kernel  at  a  smooth  point, 

/Si  Kj  (Qi,  Pj)  dSQi  =  2rc,  i  =  1,2  (97) 


to  write  (94)  and  (95)  in  the  form 

4lroi(Pl)  ■  fs, tai(Q,)  -  °,(p,)  dSQ1 

4 a2(Q2>  ^  dSQ2  =  Ul  <98) 

4TO2<P2>  -  fS2  [02(Q2)  '  °2(P2)  dS« 

4  0,(Q,)  3^  ^  dSQ1  =  U2  3^  (99) 

It  was  shown  by  Landweber  and  Macagno  [16]  that  the  singularity  of  the  ker¬ 
nel  is  removed  by  the  foregoing  procedure. 

After  the  integral  equations  have  been  solved  for  the  source  strengths 
c^i  and  a12  when  U|  =  1  and  U2  =  0,  and  for  a21  and  a22  when  U1  =  0  and  U9  = 
1,  where  Ojj  denotes  the  source  distribution  on  the  j-th  body  due  to  the 
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motion  of  the  i-th,  then,  for  the  2-D  case,  the  added  masses  can  be  obtained 
from  the  generalized  Taylor  formulas.  For  the  2-D  case,  we  obtain 

=  2nppc^  x  CTj|dSj  -  pA-j,  A22  —  2jcp ^^22d^2  ”  P^2 
where  A|  and  A2  are  the  areas  bounded  by  C|  and  c2,  and 

A|2  =  A2-j  ~  2icp^^  xa12ds2  —  2?cp|)^  xo2jds2 

and  for  the  3-D  case, 

Aji  =  4*P£sj  xandsi  ■  pVp  A.22  =  47Ip£s2  xa22ds2  '  P^2 
where  VI  and  V2  are  the  volumes  of  the  interiors  of  and  S2,  and 
Aj2  =  A2j  =  4jtp^g2  xa12dS2  =  47tp£Si  xo2|dS| 

2.  Elimination  of  peaks  of  integrands 

When  the  gap  between  the  bodies  is  small,  then,  at  points  P1Q2  and 
P2Q1  in  the  neighborhood  of  the  gap,  rpQ  will  also  be  small  and  the  integrand 
will  have  a  sharp  peak.  There  are  two  ways  to  evaluate  such  an  integral  accu¬ 
rately.  One  is  to  eliminate  or  reduce  the  magnitude  of  the  peak,  the  other  is 
to  use  a  quadrature  formula  which  concentrates  enough  points  in  the  neigh¬ 
borhood  of  the  peak  to  yield  the  desired  accuracy.  A  combination  of  these  two 
procedures  was  used  in  [2]  for  the  case  of  a  pair  of  circles,  and  the  resulting 
values  of  the  added  masses  agreed  very  well  with  those  from  the  method  of 
images;  but  the  procedure  for  modulating  the  peak  depended  upon  the  exact 
integral  of  the  kernel  for  a  circle,  which  is  not  available  with  other  shapes  and 
removed  only  the  peak  associated  with  the  points  P  at  0i  =  0  or  02  =  n  for  a 
pair  of  circles. 

A  procedure  suitable  for  arbitrary  shapes  is  based  on  the  property  of  the 
transpose  of  the  2-D  kernel  K(Pi,Q2), 

a 

rC2  K(Q2,Pi)  dsQ2  =$>C2  In  rpiQ2  dsQ2  =  0  (100) 
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since  the  integral  may  be  interpreted  as  the  flux  through  a  dosed  curve  C2  due 
to  a  unit  point  source  at  the  exterior  point  Pi.  Similarly  we  have 

JC1  K(Qi,P2)  dsQi  =  0  (101) 


and  for  the  3-D  kernels 

JSl  K(Qi/P2)  dSQi  =  JS2  K(Q2/Pi)  dSQ2  =  0  (102) 


Thus,  for  the  integral  over  ci  in  (87)  and  (92),  we  may  write 

d  d 

ai(Qi)  In  rp2Qi  dsQi  =  [ai(Qi)  In  rp2Qi  - 

'  3 

oi(Pi)  — —  In  rp2Qi]dsQi  (103) 

where  P^  is  the  point  on  ci  dosest  to  P2,  which  is  such  the  line  P^  P2  lies  along 
the  outward  unit  vector  npi/,  normal  to  ci  at  P  ;  see  Fig.  4.  Then,  when  Qi 
coinddes  with  P|,  the  last  integrand  becomes 


°i<pi)  a  a  ai 0Pj)  aaa 

rp2Pi'  ^np2  3npi^  rp2P1  -  rp2pi'  (nP2  +  nPi  ) 


glgj) 

rP2PT 


/-t  ^  ^  \ 

(1  +  npr  •  np2)  = 


Oi(Pp 

rP2PT 


(1  +  COS  \|/) 


where  np2  is  the  outward  unit  normal  vector  at  P2  on  c2  and  is  the  angle 
between  the  normals.  In  the  range  in  which  rp2pr  is  very  small,  7t  -  \\r  is  also 
small,  and  hence  1  +  cos  \\f  appears  to  be  small  of  second  order,  with  rp2pi' 
considered  as  small  of  first  order.  Hence  the  introduction  of  the  transpose  in 
(103)  reduces  the  magnitude  of  the  peak  to  first  order  of  smallness.  Similarly, 
the  second  integrals  of  (86),  (95)  and  (94)  may  be  written  as 


c2  3npi 
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(104) 


,  a 

02(P2)  In  rpiQ2]  dsQ2 

Ic  cri(Qi)  ~  — —  dSQi  =  Js  [ai(Qi)  -  — — 

S1  0np2  rP2Ql  8np2  rP2Ql 


‘5l<pi>a^-'p2Q.ldSQ1 

.  a  i  .  a  i 

Jc  CT2(Q2)  - - “ - dSQ2  =  Ic  [<J2(Q2)  “  I 

S2  0npi  rPlQ2  ^  S2  3npi  rPlQ2 


(105) 


°2(p^^lds® 

in  which  P2  is  the  point  on  02  closest  to  Pi. 


(106) 


3.  The  MAOF  with  nonuniform  intervals 

The  second  phase  of  treating  the  peak  of  the  integrand  was  to  select  a 
quadrature  formula  which  concentrates  points  in  its  neighborhood.  For  this 
purpose,  a  transformation  was  introduced  which  permitted  the  MAQF  to  be 
applied  with  nonuniform  integrals,  such  that  the  desired  concentration  of 
points  about  the  peaks  was  obtained.  This  will  now  be  presented. 

Let  F(9)  be  a  periodic  function,  with  period  2k,  which  has  a  sharp  peak 
for  small  values  of  0.  We  wish  to  evaluate  the  integral 

I  =  F(0)  d0 


If  the  MAQF  were  applied  directly,  a  very  large  number  of  abscissas  would  be 
needed  to  obtain  enough  points  over  the  range  of  the  peak,  since  uniform  in¬ 
tervals  are  required.  To  resolve  this  dilemma,  we  define  a  variable  co  by 

0  =  co -sin  co,  0<co<2k  (107) 

co 

Then  d0  =  (1  -  cos  co)  dco  =  2  sin2  y  dco 

Hence 
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I  =  2  Jgrt  G(©)  sin2 1  d©,  G(co)  =  F(0) 


(108) 


Clearly,  the  new  integrand  is  also  cyclic.  Hence  we  may  apply  the  MAQF  (the 
trapezoidal  rule)  by  using  uniform  intervals  A©.  This  concentrates  many  in¬ 
tervals  near  0  =  0  and  2k  where  A0  =  ^©2A©,  and  fewer  near  0  =  co  =  n,  where 

A0  =  2Aco,  as  was  desired.  Here  we  assumed  that  0  =  0  and  2k  at  the  points  of 
minimum  distance  on  both  bodies.  For  the  two-circle  case,  where  the  points 
of  minimum  distance  are  at  0j  =  0  and  02  =  k,  we  define  02  =  k  -  02  with  which 
the  above  formulas  may  be  directly  applied. 


For  small  values  of  co,  (107)  gives  0  =«  co3/6.  This  suggests  that  the 
concentration  of  points  near  0  =  0  and  2k  could  be  increased  by  removing  the 
co3-term.  This  can  be  accomplished  simply  by  modifying  the  0(co)  function  in 
(107)  by  subtracting  cd3/6,  but  that  would  violate  the  MAQF  condition  that  the 
integrand  be  cyclic.  Instead,  we  define  0(co)  as 


1  i 

0  =  co  -  sin  co  -  g  sin-5© 


(107a) 


Then 


d0  =  (2  sin2  ^  ^  sin2©  cos  ©)d© 


A  © 

=  2  sin4  2  (2  +  cos  ©)  d© 


and  hence 
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r  ,  © 

I  =  2  I  G(©)(2  +  cos  w)  sin4  ^  d©  (108a) 

which  also  has  a  cyclic  integrand.  At  small  values  of  ©,  we  now  have  A0= 

g  ©4A©,  so  that  the  concentration  of  0-intervals  near  0  =  0  is  greatly  increased 

when  A©  is  constant.  At  0  =  ©  =  k,  both  (107)  and  (107a)  give  d0  =  2d©.  The 
asymptotic  form  of  (107a)  is  0=3©5/4O,  which  can  also  be  eliminated  by 
subtracting  (3/40)  sin5©.  This  yields  the  third  transformation 
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0  =  co-  sino)  -  £  sin3o)  -  ^  sin5o> 


(107b) 


2k 

and  I  =  J  G(o)(8  +  9  cos  0)  +  3  cos2co)  sin6  ^  do)  (108b) 

The  procedures  for  obtaining  a  succession  of  nonuniform  intervals,  with 
which  the  MAQF  can  be  applied,  is  now  evident. 

Since  0  =  0  when  co  =  0  and  0  =  k  when  co  =  7t,  and  the  co-intervals  are 
uniform,  an  indication  of  the  point  concentration  is  given  by  the  value  0  =  0C 

when  co  =  n/2,  i.e.  at  half  the  number  of  points  between  0  and  n.  Substitution 
of  co  =  n/2  in  (107)  and  (107a)  gives  0C  =  33°,  23°,  and  19°  respectively.  Also,  at 
these  values  of  co  and  0C  we  have  d0  =  dco,  indicating  that  the  increments  A0 
change  from  concentration  to  dispersion  at  0  =  0C. 

With  any  member  of  this  family  of  modified  MAQF  transformations, 
the  integrand  vanishes  at  0  =  0  and  2;t,  thus  eliminating  the  peak  of  the 
original  integrand  F(0)  at  that  point.  This  property  also  reduces  the 
amplitude  of  the  peaks  occurring  at  small,  nonzero  values  of  0,  but,  as  will  be 
seen,  not  sufficiently  so  that  the  transpose  correction  for  eliminating  peaks  of 
the  kernel  would  become  unnecessary. 

When  co  is  small,  significant  figures  may  be  lost  in  computing  0  directly 
from  (107, 107a,  107b).  For  (107),  this  can  be  avoided  by  writing 

co3  co5  co7  co3  co2  co4  co6 
9  =  3!  "  5!  +  7!  "  =  6  ( 1  *  20  +  840 ' 60,480)  +  e(to) 

where  I  e(co)  I  is  less  than  the  next  term  of  the  series, 

lei  <  co  11/11!  <  1.6  xlO'7©11 

and  co3/ 6  is  computed  in  scientific  notation.  For  (107a)  and  (107b),  we  apply 
1  1 

sin3co  =  4  (3sinco  -  sin3co),  sin5co  =  (10  sinco  -  5  sin3co  +  sin5co) 

and  the  Taylor  series  of  sines  of  co,  3co  and  5co  to  obtain,  respectively. 
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(107c) 


0  = 


3co5 

40 


1-15  2”  (-l)n 


co2n“4 

(2n+l)! 


(32n'2-l) 


0  = 


5co7 


112 


105 

8 


2^  (-1)" 


CO2  n_6 
(2n+l)! 


(2-32n_1  +  52n‘2) 


(107d) 


4.  Evaluation  of  derivatives  of  added  mass 

Both  added  masses  and  their  derivatives  with  respect  to  the  gap  g  be¬ 
tween  the  cylinders  occur  in  the  expressions  for  the  forces  on  the  bodies.  In 
Ref.  [1],  these  derivatives  were  obtained  by  numerical  differentiation  of  the 
smoothed  added-mass  data,  with  unsatisfactory  results  at  small  gaps,  because 
the  derivatives  were  approaching  infinity  and  small  errors  in  the  data  were 
amplified. 


The  accuracy  of  the  results  from  numerical  differentiation  of  the 
added-mass  coefficients  has  been  greatly  improved  by  a  simple  and  important 
change.  Since  the  error  depends  on  the  magnitude  of  the  second  derivative, 
this  suggests  that  the  asymptotic  law  of  the  added-mass  be  applied  to  reduce 
the  error.  For  example,  according  to  (66),  ki(g)  varies  linearly  with  Vg  (or  the 
parameter  £  or  y  of  (62))  at  small  values  of  g,  as  is  shown  in  Fig.  2.  Then 
dkj/dy  can  be  obtained  by  numerical  differentiation  with  greatly  reduced 
error,  with  which,  by  (62),  we  get 


dki  _  dki  dy  _  1  /2s  >  1/2  dkt 
dg  "  dy  dg~2[p2gJ  ^ 


(109) 


without  further  loss  of  accuracy. 


The  parameter  y  is  also  useful  for  improving  the  accuracy  of  the  added- 
mass  derivatives  for  spheres.  This  is  not  indicated  directly  by  the  asymptotic 
formula  (78),  since  the  rate  of  variation  of  the  added  masses  at  very  small 
values  of  £  or  y  is  much  greater  than  at  its  moderately  small  values 
corresponding  to  the  practical  range  for  g/a  >  0.01  (for  equal  spheres)  due  to 
the  factor  In  y.  Since  dk|/dy  varies  slowly  in  that  range,  it  could  be  evaluated 
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with  sufficient  accuracy  from  the  computed  values  of  the  added  masses. 
Their  derivatives  with  respect  to  g  were  then  obtained  from  (109). 

To  obtain  the  numerical  results  for  the  derivatives  (other  than  the 
'exact'  ones)  shown  in  the  following  tables,  Lagrange's  five-point 
interpolation  formula  was  used,  without  smoothing  the  calculated  values  of 
the  added-mass  coefficients. 

5.  Numerical  results  for  pairs  of  circles 

a)  Exact  results 

All  the  numerical  results  were  obtained  on  an  IBM  RISC/6000 
minicomputer.  First,  an  accurate  set  of  added  masses  for  a  pair  of  equal  circles 
was  obtained  by  summing  infinite  series  of  doublet  strengths.  These  were 
calculated  by  means  of  the  basic  recurrence  formulas,  such  as  (3)  and  (6),  and 
by  the  various  forms  of  the  Herman  formulas.  Of  these,  the  parametric  forms 
in  (39b,  40b,  41b)  were  found  to  be  most  efficient  and  their  efficiency  was 
greatly  enhanced  by  applying  the  truncation  corrections  given  in  (49,  53,  61)  in 
summing  the  series  for  the  added-mass  coefficients.  This  is  demonstrated  in 
Table  2  by  calculations  of  k\2  for  equal  circles,  where  R12  is  given  by  the 
trapezoidal  integral  (49).  This  table  supplements  Table  1  which  corresponds 
to  the  zero  gap.  We  see  that  the  approximate  ratios  of  the  number  of  terms  N 
which  give  the  same  accuracy  without  and  with  the  truncation  correction  are 
10  at  g/a  =  0.0001,  4  at  0.001,  and  2.5  at  0.01.  This  shows  that  the  truncation 
correction  is  very  effective  at  very  small  gaps,  and  that  its  effectiveness 
diminishes  with  increasing  gaps.  By  (62),  the  parameter  C,  is  related  to  these 
gaps  by  £  =  2^ g/a  =  0.02,  .0632,  0.20,  respectively.  Thus  £  is  not  very  small  at 
g/a  =  0.01,  and  the  series  for  ki2  at  larger  gap's  (say  g/a  >  0.10)  converge  fast 
enough  so  that  little  is  gained  by  applying  the  truncation  correction. 

It  was  found  that,  for  g/a  >  0.01,  the  added-mass  coefficients  obtained 
from  the  recurrence  and  the  various  forms  of  the  Herman  formulas  were 
identical  to  the  sixth  decimal  place.  This  would  not  be  the  case  at  much 
smaller  gaps,  corresponding  to  very  small  values  of  £,  when  N  is  very  large, 
because  of  accumulation  of  round-off  errors  in  using  the  recurrence  formulas. 
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For  applications  requiring  accurate  calculations  at  very  small  gaps,  we  observe 
that 


sxnh 

%  = — 7~rz  =  2[cosh  (n-l)C  +  cosh  (n-3)C  +  ...  +  cosh  £],  n  even 
sinh£ 

tn  =  2[cosh  (n-l)£  +  cosh  (n-3)£  +  ...  +  cosh  2Q  +  1,  n  odd  (110) 

which  can  be  readily  derived  by  expressing  the  hyperbolic  -functions  in  their 
exponential  forms.  This  would  be  useful  when  £  is  so  small  that  tn  could  not 
be  calculated  accurately  from  its  definition,  even  with  double  precision. 

Values  of  ki,  ki2  for  g/a  >  0.01,  designated  as  'exact',  will  be  given  in 
Table  4.  'Exact'  values  of  ki  and  k2  are  not  listed  separately  in  the  table 
because  these  were  identical  to  the  calculated  values  to  the  six  decimals 
shown.  The  columns  for  ki  and  ki2  are  actually  bkj  and  bki2,  but  it  is 
convenient  to  take  b=l  hereafter.  The  'exact'  added-mass  derivatives  given  in 
this  table  were  also  computed  both  by  the  recurrence  formulas  (18)  and  by  the 
derivatives  of  the  Herman  formulas  given  in  (67,  68,  69).  These  agreed  to  the 
six  decimals  shown  in  Table  4. 

It  is  also  of  interest  to  show  the  variation  of  ki  and  ki2  with  £.  This  is 
depicted  in  Fig.  3  in  which  the  asymptotic  linearity  with  £  or  y,  according  to 
(52,  58,  66),  is  clearly  shown. 

b)  Computed  results 

For  the  case  of  two  circles  of  radii  a  and  b,  the  kernels  of  (88)  are 
constant  and,  according  to  (89),  the  constants  are  given  by 

K(Pi,Qi)=^  K(P2,Q2)=^ 


The  first  integrals  in  (86)  and  (87)  are  then  proportional  to  the  total  source 
strength  on  each  circle,  which  is  zero  by  the  Gauss  flux  theorem.  The  kernel 
of  the  second  integral  of  (86)  is  given  by 


K(Pi,Q2)=r^-in 

dnpi 


rPlQ2  = 


1  drpiQ2 
rPlQ2  dnpi 
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Here 


rpiQ2  =  [a2  +  b2  +  c2  +  2bc  cos  0Q2  -  2ca  cos  9pi  -  2ab  cos  (0Q2-0pi)]1//2 

and  3rpiQ2/3npi,  the  cosine  of  the  angle  between  the  vector  from  Q2  to  Pi  and 
ni,  can  be  obtained  by  differentiating  rpiQ2  with  respect  to  a.  This  gives 

K(PlrQ2)  =  ~2 - [a  -  c  cos  0pi  -  b  cos  (0Q2  -  0pi)]  (111) 

rPlQ2 


Similarly,  we  obtain 


K(P2,Qi)  =  - - in  rp2Qi  = 

dnp2 


[b  +  c  cos  0p2  -  a  cos  (0p2  -  0qi)1  (112) 


rP2Ql 


and  for  their  transposes. 


d  1 

K(Q2,Pi)  =  — —  in  rpiQ2  =  ~y~ 


3nQ2 


rPlQ2 


a  1 

K(Qi,P2)  =  - - in  rp2Qi  =  y 

3nQi  ^  / 


P2Q1 


[b  +  c  cos  0Q2  -  a  cos  (0pi  -  0q2)]  (113) 


[a  -  c  cos  0p2  -  b  cos  (0qi  -  0p2)3  (1 14) 


The  integral  equations  (86)  and  (87),  with  the  transposes  of  the  kernels 
applied  to  eliminate  the  peaks  as  shown  in  (103)  and  (104),  now  become 

tcoi(Pi)  +  j  [a2(Q2)  K(Pi,Q2)  -  c(PJ  K(Q2,Pi)]  bd02  =  Ui  y-  (115) 

c2  2  0nPl 

'  dx 

7to2(P2)  +  J  [cji(Qi)  K(P2,Qi)  -  o(P  )  K(Qi,P2)]  ad0i  =  U2  r -  (116) 

C1  1  »nP2 

/ 

Here  the  point  P^  is  such  that  rp2pp  is  the  minimum  distance  from  the  point 

/ 

P2  on  one  body  to  the  other;  see  Fig.  4.  For  this  case,  P^  is  the  point  of 

intersection  of  the  line  from  P2  to  the  center  of  the  other  circle,  and  one  can 
readily  show  that  0pi'  and  0p2'  are  given  by 


41 


(11 7) 


tan  0pi'  = 


b  sin  0p2 
c  +  b  cos  0P2 


tan  0p2'  = 


a  sin  0pi 
c  -  a  cos  0pi 


The  integrals  were  discretized  by  applying  the  modified  MAQF  of  (107) 
and  (108)  to  obtain  a  set  of  2N  linear  equations  for  the  o's.  For  integration 
over  circle  A,  in  which  the  peak  occur  near  k  instead  of  zero,  02  in  (107)  must 

9  9  9 

be  replaced  by  its  supplement  0^  i.e.,  =  jc  -  02  =  to  -  sin  co.  Then  0^  and  co  vary 

from  0  to  2n  as  02  varies  from  k  to  -k.  The  discretized  forms  of  (115)  and 
(116)  then  become,  in  the  form  of  an  iteration  formula. 


(n+l)  TT 
7iaH  =  Ui 


dx  47tb 
3nii  N 


IN[a2j 

j=l 


(n) 


(n) 


5S 


Kii2j  -  a2.,  K2jii]  sin2 


(n+l) 

2i 


U23n2l  N  fa  1<Jij 


(n) 


2«i 


K2iij  -  Kij2i]  sin2  -g 


(118) 


where  the  superscripts  (n)  and  (n+l)  indicate  the  nth  and  the  (n+l)th 

9  / 

iterations.  The  values  on'  and  G2i',  corresponding  to  a(PJ)  and  a(P2),  were 
obtained  by  linear  interpolation  between  the  values  of  a  at  the  extremities  of 

9  9 

the  co-interval  in  which  P  or  P2  lies.  The  initial  approximations  to  start  the 
iteration  were  taken  to  be 


(0)  1  _  3x 

ali  =_U1  ^ - 

11  n  3nii 


(o)  i  TT  ax 

°2i  =~U2^ — 
^  n  an2i 


(119) 


/ 

Graphs  of  the  source  strengths  ai(0i)  and  02(0^  and  of  the  original 

integrand  02(82)  K(Pi,Q2)  of  (115)  are  given  in  Figs.  5  and  6  for  the  case  Ui  =  1, 
U2  =  0,  and  a  =  b.  Results  for  three  gaps,  g/b  =  0.01,  0.05,  0.10  in  Fig.  5  show 

9 

that  both  oi  and  02  peak  sharply  at  01  =  0^  =  0,  with  peak  maxima  of  01  (0)  =  1.7 


and  02(0)  =  1.5  at  g/b  =  0.01  and  the  smaller  values  Oi(0)  =  0.80,  02(0)  =  0.64  at 
g/b  =  0.05,  and  oi(0)  =  0.60,  02(0)  =  0.44  at  0.10,  although  the  peaks  of  02  are 
consistently  sharper  than  those  of  gj.  As  the  gap  approaches  infinity,  the 
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1 

source  distributions  become  Oi  »  -  cos  0,  oi(0)  «  l/7t  and  02  =  0.  At  0i  =  0  =  tc, 

n  2 

the  limiting  source  strengths  are  oi(tc)  =  -  1/rc  and  02(71)  =  0,  and  the  values  at 
the  small  gaps  are  negative  but  deviate  little  from  the  limiting  values  at  g  =  «>. 

Figures  6a,  6b,  6c  show  curves  for  g/b  =  0.01  corresponding  to  three 
fixed  points  0pi  which  give  0'p2'  =  0,  6.5, 11.9  degrees.  We  see  that  the  peaks  of 
the  original  integrands  are  an  order  of  magnitude  higher  than  those  of  the 
source  strengths  at  this  gap.  For  example,  the  maximum  value  of  the  kernel 

K(Pi,Q2)  at  0i  =  0  occurs  at  0^  =  0  and  has  the  value  1/g  =  100  versus  the  peak 

02(0)  =  1.5,  in  good  agreement  with  the  calculated  peak  value,  02(0)  K(pi,Q2)  = 
151.6.  At  0'p2'  =  6.5  and  11.9  degrees,  the  peak  values  of  the  original  integrand 
are  27.1  and  4.0,  much  less  than  the  magnitude  of  the  central  peak. 

With  the  transpose  correction  and  the  associated  source  strength 

^2(0 2p'),  Pea^s  were  essentially  eliminated.  When  02(71)  was  used  instead 

/ 

of  O2(02p/),  as  is  done  in  Ref.  [2]  for  all  points  Pi,  the  peaks  were  overcorrected 

and  the  peaks  of  greater  magnitude  but  of  opposite  sign  shown  in  Figs.  6b  and 

6c  were  obtained.  This  is  a  direct  consequence  of  the  sharp  peak  of  the  02- 

curve  in  Fig.  5a,  which  shows  that  02(0)  is  more  than  twice  the  values  of  02  at 
/ 

02p,  =  6.5  and  11.9  degrees.  This  shows  that,  although  the  kernel  peaks  are 

much  larger  than  those  of  the  source  distributions,  the  latter  play  an 
important  role  in  the  elimination  of  the  peaks  of  the  integrand. 

In  Fig.  7,  the  ordinate  scales  of  Fig.  6  are  magnified  to  show  the  fine 
structure  of  the  nearly  flat  curves  of  the  integrand  when  its  peaks  have  been 

t 

removed.  The  initial  dip  of  the  curve  for  02p,  =  0  is  due  to  a  combination  of 

the  small  value  of  the  gap,  g/b  =  0.01,  and  the  associated  very  sharp  peak  of 
/ 

02(82);  this  was  verified  analytically  by  applying  the  expressions  for  the  kernel 
in  (111)  and  its  transpose  in  (113).  For  the  other  two  values  of  02p/,  many 
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more  points  would  have  been  needed  to  show  smooth  curves  of  their 
oscillations;  but  since  these  curves  were  not  used  quantitatively,  and  the 
purpose  of  the  figure  was  to  show  the  need  to  concentrate  many  points  of  a 
quadrature  formula  in  the  gap  region  to  take  this  fine  structure  into  account, 
the  present  figure  sufficed.  This  was  accomplished  by  using  the  modified 
MAQF  (107)  to  concentrate  points  in  the  regions  of  the  peaks  and  increasing 
the  number  of  points  as  the  gap  decreases.  In  Ref.  [2],  three  Gauss  quadrature 
formulas  of  order  20,  two  of  which  extended  over  the  peak  region,  were  used 
for  the  same  purpose. 

The  foregoing  hypotheses  were  tested  by  running  the  numerical 
experiments  of  Table  3.  Condition  1  gave  'exact'  values;  comparison  of 
conditions  1  and  4  confirms  the  importance  of  removing  integrand  peaks, 
although  the  accuracy  of  4  improved  greatly  at  N  =  100.  Condition  2 
simulates  the  above-described  technique  of  Ref.  [2],  in  which  oi(0)  and  02(0) 
are  used  instead  of  ai(0ipO  and  O2(02pO.  The  accuracy  improved  greatly  by 
increasing  N  from  40  to  100  and  the  latter  result  agrees  very  well  with  that  of 
Ref.  [2]  given  in  Table  5.  Condition  5  is  the  procedure  used  in  Ref.  [1],  which 
required  N  =  600  to  obtain  acceptable  accuracy.  Thus  it  appears  to  be  most 
efficient  and  accurate  to  remove  all  the  peaks  of  the  integrand  and  to  use 
some  means  of  concentrating  points  in  the  peak  regions. 

The  present  results  are  compared  with  the  'exact'  ones  in  Table  4  and 
with  those  of  Ref.  [2]  in  Table  5.  The  latter  shows  that  the  values  of  kj  and  ki2 
given  by  the  present  methods  agree  with  the  exact  values  to  six  decimals.  In 
Ref.  [2],  the  results  were  correct  to  only  two  decimals  at  gaps  of  0.01  and  0.02, 
three  decimals  at  0.03  and  0.04,  and  to  at  least  four  decimals  at  0.05  and 
beyond.  The  greater  accuracy  with  the  present  methods  is  evident. 

The  need  for  high  accuracy  becomes  apparent  when  the  data  are  used  to 
obtain  the  derivatives  of  the  added  masses  by  numerical  differentiation, 
which  depend  upon  differences  between  successive  values.  From  Table  5,  we 
see  that  the  differences  between  values  at  gaps  of  0.01  and  0.02  for  ki,  0.0313 
from  Ref.  [2],  versus  the  exact  value  of  0.0337,  is  in  error  by  7  percent.  Clearly, 
then,  numerical  differentiation  could  not  be  used  in  Ref.  [2].  Instead,  the 
derivatives  of  the  original  integral  equations  were  solved  for  the  derivatives 
of  the  source  strengths,  which  then  yielded  the  derivatives  of  the  added- 
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masses  directly.  This  procedure  gave  the  values  shown  in  Table  6,  which 
have  an  error  of  2.5  percent  at  g/a  =  0.01,  but  are  very  accurate  of  gaps  of  0.02 
and  beyond.  Since  the  peaks  of  the  integrands  of  the  differentiated  integral 
equations  were  much  higher,  the  order  of  the  two  Gauss  quadrature  formulas 
in  the  symmetrical  peak  regions  was  increased  from  20  to  40.  The  values  of  kj 
and  ki2,  obtained  with  the  four  Gauss  quadrature  formulas,  each  of  order  40, 
are  not  listed  in  Ref.  [2].  The  present  method  gives  more  accurate  results  at 
g/b  =  0.01  and  0.02,  but  at  0.03  and  beyond,  both  are  highly  accurate.  The 
present  method,  however,  attains  the  accuracy  with  a  much  simpler 
computer  program,  requiring  the  simultaneous  solution  of  two,  instead  of 
four,  integral  equations,  with  the  kernels  of  the  two  additional  equations 
having  peak  magnitudes  of  order  l/g2  versus  1/g  for  the  original  pair.  The 
numerical  results  by  the  present  method  were  obtained  with  a  minicomputer; 
those  of  Ref.  [2]  for  ki  and  ki2  required  a  supercomputer. 

Results  for  unequal  circles  are  given  in  Tables  7,  8  and  9  for  a/b  =  4, 
16,«>o.  'Exact'  results  are  not  shown  separately  since  these  agreed  exactly  with 

the  computed  results  to  the  number  of  digits  shown.  This  was  accomplished 

a  a 

by  using  the  modified  MAQF  (107)  on  both  circles  and  taking  N(^)  =  ^  N(l). 

Actually,  it  was  necessary  to  increase  the  number  of  points  only  on  ci,  and 
there  only  in  the  peak  regions,  as  we  have  verified  by  using  (107b)  to  increase 
the  point  concentration  there  for  a/b  =  16,  with  a  smaller  value  of  Nj.  The 
results  in  Table  7,  however,  were  obtained  with  N  =  120  on  both  circles,  and 
those  in  Table  8  with  N  =  600  on  both  at  g/b  =  0.01,  although  40  points  on  the 
smaller  circle  would  have  matched  the  accuracy  attained  with  a/b  =  1.  Fewer 
points  were  required  at  larger  gaps. 

The  coefficients  ki  in  Tables  7  and  8  are  seen  to  be  given  approximately 
by  (a/b)2  followed  by  six  decimals.  As  is  seen  from  (39b),  however,  (a/b)2  is  a 
known  additive  constant,  and  the  remaining  terms  need  to  be  calculated 
accurately  to  only  the  usual  number  of  seven  significant  figures.  Similarly,  in 
solving  the  integral  equations  for  the  source  distribution  on  the  larger  circle, 
the  known  contribution  to  the  distribution  from  the  doublet  of  strength  a2  at 
the  center  of  the  large  circle,  a0  =  (1/rc)  cos  0i,  would  also  subtract  out  the 
additive  constant  (a/b)2  from  the  values  of  ki  and  require  only  the  same 
seven-digit  accuracy  for  the  remainder.  This  procedure  can  be  generalized  to 


45 


apply  to  a  noncircular  cylinder  by  obtaining  first  the  source  distribution  aQ  on 
the  surface  of  that  cylinder  when  the  other  cylinder  is  not  present,  and  then 
solving  the  pair  of  integral  equations  for  the  interaction  source  distribution  a- 
aQ. 


Graphs  of  the  force-coefficient  data  in  Tables  4,  7,  8  and  9  are  shown  in 
Figs.  8a  to  8d.  Here  Cfi  denotes  the  coefficient  of  the  force  on  the  circle  of 
radius  a  when  that  circle  has  the  velocity  Ui  and  the  circle  of  radius  b  is  at 
rest,  and  Cf2  that  for  the  force  on  the  circle  of  radius  b.  These  coefficients  are 
defined  by 

Cfi  =  "  Fl  ~2  Cp2  = 

pb  U| 


F2 

pbU? 


where  the  subscript  T  corresponds  to  the  a-circle  and  '2'  to  the  b-circle. 

Similarly,  when  the  a-drcle  is  at  rest  and  the  b-circle  is  moving  with  velocity 

*  * 

U2,  the  forces  are  designated  as  Fj  and  F2  and  the  coefficients  are  defined  by 


For  these  special  cases,  we  obtain  from  (15)  and  (16)  with  b  =  1 


1  .  1  .  *  1  . 

Cfi  =  2  K  kj,  Cf2  =  ft  (^  ki  +  ki2),  Cp1  =  -  n  (ki2  +  2  k2),  Cp2  =  -  rck2  (120) 


Although  the  procedures  for  computing  the  results  for  a  =  °°  are  not 

presented  until  the  next  section,  it  is  convenient  here  to  consider  this  case  as 

the  limit  of  the  sequence  of  circles  of  increasing  radii.  Curves  of  Cfi  and  Cf2 

do  not  appear  in  Figs.  8a  and  8b  since  motion  of  the  wall  towards  a  stationary 

cylinder  would  require  an  infinite  added  mass.  We  see  that  the  forces  are 

always  repulsive  and  approach  infinity  as  the  gap  approaches  zero.  Values  of 
*  * 

Cp-j  and  Cp2  are  not  shown  explicitly  in  Table  4  for  a/b  =  1  since,  by  symmetry. 
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*  *  * 
we  have  Cfl  =  -  Cf2  and  Cp2  =  -  Cfi-  The  graphs  show  that  Cfi  and  CF2  vary 

* 

monotonically  with  a/b  for  all  values  of  g/b.  Cp2  and  CF1  also  vary 
monotonically,  but  only  for  values  of  g/b  from  0  to  0.16  and  0.04  respectively. 

When  a/b  =  1,  Figs.  8a  and  8b  and  Table  4  show  that  the  force  on  the 
stationary  cylinder  is  about  five  times  larger  than  that  on  the  moving  one  at 
g/b  =  1,  and  about  50  percent  larger  at  g/b  =  0.10.  This  indicates  that  it  would 
be  easier  to  detect  hydrodynamic  interaction  effects  experimentally  on  the 
stationary  cylinder,  for  which  the  predicted  value  Cf2  =  0.75  at  g/b  =  0.50  is 
large  enough  to  be  readily  measured.  For  other  diameter  ratios  at  g/b  =  0.50, 
the  force  coefficient  on  the  large  cylinder  is  greater  no  matter  which  of  the 
cylinders  is  moving. 

6.  Results  for  a  circular  cylinder  (or  a  sphere)  and  a  wall. 

The  case  of  a  circular  cylinder  moving  along  the  normal  to  an 
impervious  wall  is  the  limiting  case  of  a  pair  of  circles  as  the  radius  of  one  of 
the  circles  approaches  infinity.  Each  of  these  large  circles  is  taken  to  be  at  rest 
as  a  so  that  g/a  -»  0,  in  contrast  to  the  results  for  a/b  =  16  in  Table  8, 

which  are  given  to  g/a  =  6.25.  Thus  results  for  the  present  case,  given  in  Table 
9,  should  resemble  those  of  Table  8  only  for  small  gaps,  say  g/b  <  1.0,  at  which 
g/a  <1/16. 

As  is  well  known,  an  alternative  approach  is  to  replace  the  wall  by  the 
mirror  image  of  the  circle  in  it.  This  method  transforms  the  problem  to  the 
case  of  a  pair  of  equal  circles  for  which  results  are  given  in  Table  4.  Both 
approaches  will  now  be  treated. 

As  a  approaches  infinity  in  (40b)  and  (41b),  we  obtain  in  the  limit 

k!=1+2I"«=»-1'  k2=2i2  (121) 

ki2  =  -  2X2(9,  kl2  =  -k2  (122) 
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Then,  by  (15)  and  (16), 


Fi  =  -  (A12  +  \  A22)  =  -  7cpb2(k12  +  \  k2)U2  =  7cpb2i2U2  (123) 

F2  =  - 1 A 22  U2  =  _  \  Jtpb 2k2  =  -  7cpb2t2U2  (124) 


With  b  =  1,  these  give  the  force  coefficients 


Cfi  =  — ^7  =  Cf2  =  -  tcZ2(0 

pbU2 


(125) 


The  value  of  £  for  given  values  of  a,  b,  and  c  can  be  obtained  from  (63), 
which  may  be  written  as 


cosh;  =  i+|+^f+^ 


In  the  present  case  of  a  =  «>,  we  get 

cosh  £-1=2  sinh2  2  =  b  C  =  2  sinh'1  ^bj^2 

To  obtain  Z2,  we  first  apply  (69)  for  a  ->  «,  to  get 

1O0  d  -i 

d£Xn 


v '  -  V°°A  T‘2  _  .  2  y°°  t'3  r' 
l2  ArX n  “  2  Tn  xn 


As  a  and  c  approach  infinity  in  (69),  this  gives 

x2©=^(0/c-=-5cSchcs7t;3t; 


(126) 


The  force  coefficients  (125),  for  the  case  of  a  cylinder  moving  normal  to  a  wall 
then  become 
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(127) 


cschC^t;3*; 

Next  consider  the  two-equal-circle  approach.  By  (14),  the  kinetic  energy 
of  the  fluid,  with  Ui  =  -  U2  and  An  =  A22,  is 

Tc  =  [A22c(gc)  *  A]2c(gc)]U2/  gc  =  2g 

and  that  for  the  half  space  with  Ui  =  0  is  given  by  2T  =  A22(g)  Ui-  Here  the 

subscript  c  indicates  values  for  the  pair  of  equal  circles,  and  the  other 
quantities  refer  to  the  half  space  when  the  wall  is  present.  Since  Tc  =  2T,  we 
obtain  the  simple  relation  A22  =  A22C  •  Ai2c  or,  in  nondimensional  form, 

k2(|)  =  k2c(f)-k12c(f)  (128) 

Hence,  differentiating  with  respect  to  g,  we  obtain 

k2(|)  =  2[k2c(f)-ki2c(f)]  (129) 

The  derivation  of  the  results  Fi  =  -  F2  in  (123)  and  (124)  was  based  on 
the  formulas  for  ki  and  k^  for  a  circular  cylinder  in  (121)  and  (122).  We  shall 
now  show  that  Ft,  the  force  on  the  wall,  is  also  given  by  Fi  =  -  F2  for  bodies  of 
arbitrary  shape. 

According  to  the  Lagally  theorem  [5],  the  force  on  a  nonrotating  body  in 
an  irrotational  flow  can  be  expressed  in  terms  of  the  strength  and  location  of 
the  singularities  within  it,  the  velocities  due  to  external  mechanisms  at  their 
locations,  and  the  acceleration  of  the  displaced  fluid  (which  is  zero  in  the 
present  case).  Since  the  image  body  and  the  wall  (considered  as  a  circle  of 
infinite  radius)  contain  the  same  set  of  singularities,  the  force  must  be  the 
same  on  both.  But,  by  symmetry,  the  force  on  the  image  body  is  equal  and 
opposite  to  that  on  the  actual  one.  Hence  that  is  also  true  for  the  force  on  the 
wall,  as  we  wished  to  show. 

Also  valid  for  bodies  of  arbitrary  shape  are  equations  (123, 124, 128, 129), 
the  derivation  of  which  did  not  depend  on  the  body  shape.  Since  Fi  +  F2  =  0, 
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we  now  find,  from  (123)  and  (124),  that  the  relation  A22  =  -  A12,  given  in  (129) 
for  a  circular  cylinder,  also  applies  to  arbitrary  shapes.  Thus,  the 
determination  of  the  forces  Fi  and  F2  on  a  body,  approaching  a  wall  at 
constant  speed  requires  only  that  the  added  masses  be  determined  by  the 
integral-equation  procedure  that  has  already  been  described,  applied  to  the 
body  and  its  mirror  image. 

For  the  case  of  a  sphere  moving  away  from  a  wall,  we  treat  the  wall  as  a 
sphere  of  infinite  radius.  Then,  as  a  approaches  infinity  in  (75)  and  (76),  we 
obtain  in  the  limit 


k2  -  2 


1  +  3Z? 


sinh3^ 


sinh3(n+l)£j 


-  2  ^3^0 ' 1 


(130) 


(131) 


Then,  by  (15)  and  (16), 


Fi=- 


A12  +  2  A22 


^  2  4  -f- 

U2  =  -3  Ttpb-5 


k12  +  2  k2 


=  7cpb3i3(C)U2 


F2  =  -  ^  A 22  U2  =  - 1  rcpb3  k2  U2  =  -  7cpb3Z3(QU2 


These  give  the  force  coefficients 

CF1  =  2Fl-2  =  =  -  Cp2  (132) 

7cpb2U2 
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As  for  the  case  of  a  circle  and  a  wall,  X3(Q  can  be  expressed  in  terms  of  k12  and 
k12.  From  (76),  we  obtain  X3(Q  =  xn ,  £3  =  -  3Z~  xn  x'  and  then 

i3  =  L|/c'  =  -|cschC-S“%  t'  (133) 

where  c'  and  x^  are  given  in  (69)  and  (70). 

7.  Ellipse-circle  combinations 

The  procedure  developed  for  pairs  of  circular  cylinders  will  now  be 
applied  to  the  case  of  an  elliptical  cylinder  approaching  a  circular  one  along 
their  line  of  centers.  Since  an  exact  solution  is  not  available  for  this  case,  we 
proceed  immediately  to  the  integral-equation  for  the  source  distributions 
ai(Qi)  and  02(Q2)  formulated  in  (86)  and  (87),  in  which  the  first  integral  of  (87) 
over  the  circle  C2  vanishes,  as  was  shown  in  Section  5. 

The  peaks  of  the  integrands  when  the  gap  is  small  were  eliminated  by 
introducing  the  transposes  of  the  kernels  and  locating  the  points  P1  and  P2 

defined  in  (103)  and  (104).  Since  C2  is  a  circle,  the  point  P2  is  again  given  by 

(117);  but  to  find  Pj,  one  must  determine  the  point  on  the  ellipse  at  which  the 

distance  from  an  exterior  point  P2  on  the  circle  is  a  minimum.  This  can  be 
found  from  the  condition  that  the  line  normal  to  the  ellipse  at  P^  passes 

through  the  point  P2(x2,y2)-  In  terms  of  the  parametric  form  of  the  equation 
of  the  ellipse,  x  =  a\  cos  y,y  =  a.2  sin\|/,  this  gives 

'  2  2  ' 
tan  \j/(P1)  =  y2/[aiX2  -  (a1  -  a2)  cos  \\f  (P-,)] 


Since  \|/  is  small  in  the  region  of  the  peaks,  this  equation  can  readily  be  solved 
by  iteration,  with  \\f  =  0  in  the  right  member  as  the  first  approximation. 

Results  for  an  ellipse  approaching  a  circle  are  given  in  Tables  10  and  11. 
In  the  former  a\  =  2b  and  a2  =  b,  in  the  latter  a\  =  b,  a2  =  2b,  where  b  is  the 
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radius  of  the  circle,  in  terms  of  which  the  added  masses  and  forces  are 
nondimensionalized.  Comparison  of  Tables  4  and  10  shows  that,  although 
the  transverse  dimensions  and  the  asymptotic  added  masses  as  g  — » °°  are  the 
same,  the  added  mass  and  forces  are  appreciably  less  for  the  ellipse-circle  pair. 
On  the  other  hand  when  Tables  4  and  11  are  compared,  it  is  seen  that  the 
interaction  effects  are  considerably  greater  for  the  ellipse-circle  pair.  These 
results  are  reasonable  since,  with  ai/a2  =  2,  the  radius  of  curvature  at  the 
leading  edge  of  the  ellipse  is  b/2,  and  with  a.\/&2  -  1/2,  is  4b.  The  accuracy  of 
the  results  in  Tables  10  and  11  was  judged  by  the  effects  of  increasing  the 
number  of  points  on  each  body.  To  obtain  the  indicated  accuracy,  N  =  300  was 
used  on  both  bodies  in  both  cases. 

8.  Results  for  two  spheres  bv  method  of  integral  equations 

Two  spheres  of  radii  a  and  b  are  moving  with  velocities  Ui  and  U2 
along  their  line  of  centers.  We  shall  use  two  spherical  coordinate  systems, 
(Rl,  0i,  Yi)  and  (R2,  02/  Y2)/  with  origins  at  the  centers  of  the  spheres  at  x  =  xi 
and  x  =  X2  for  the  spheres  of  radii  a  and  b  respectively.  With  the  x-axis  as  the 
polar  axis,  the  spherical  coordinates  are  related  to  a  fixed  rectangular 
coordinate  system  by 


x  -  Xi  =  Rj  cos  0i,  y  =  Ri  sin  0j  cos  Yi/  z  =  Rj  sin  0j  sin  Yi/  i  =  1/2 


with  X2  -  xi  =  c. 

Relative  to  the  center  of  one  of  the  spheres  of  radius  a,  we  have 

rPQ  =  RP  +  RQ-2  RPRqG(P,Q) 


where  G(P,Q)  =  cos  0p  cos  0q  +  sin  0p  sin  0q  cos  (yp~Yq) 


we  see  that  G(P,Q)  =  G(Q,P).  Then  the  kernels  of  the  integral  equation  (98) 
become 


K(P,Q)  =  -  [ 


— - — 1 

3Rp  rPQjRP  =  a 


-j—  [1  -  G(P,Q)]  = 
rPQ 


[1  -  G(P,Q)]~1/2 
2V2  a2 


(134) 
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We  see  that  the  kernel  is  symmetric;  i.e.  K  (P,Q)  =  K(Q,P). 

Next,  when  the  points  P  and  Q  are  on  different  spheres  of  radii  a  and  b, 
we  obtain 

rPlQ2  =  (c  "  xpi  +  XQ2)2  +  (ypi  -  yQ2)2  +  (ZP1  -  ZQ2)2 
=  a2  +  b2  +  c2  -  2c  (a  cos  0pi  +  b  cos  0q2)  -  2abG  (Pi,Q2)/  0'  =  ft  -  0 

where 

G(Pi,Q2)  =  -  cos  0pi  cos  0q2  +  sin  0pi  sin  0q2  cos  (\|/pi  -  \j/Q2) 

K(Pi,Q2)  =  -  ~  [a  -  C  cos  0p,  -  bG(Pi,Q2)]  (135) 

K(Q2,pi)  -  -  ~  ^  [b  -  c  cos  eQ2  -  aG(P,,Q2>]  (136) 

and  similarly 

rP2Qi  =  a2+b2+c2  -  2c(a  cos  0qi  +  b  cos  0j,2)  -  2abG  (P2/Q1) 

G(P2,Qi)  =  -  cos  0p2  cos  0qi  +  sin  0P2  sin  0qi  cos  (yp2  -  vqi) 

k(p2/Qi>  =  -  ^  ^  “  ^3^  lb  -  c  “s  4  -  aG  (p2,Ql)l  (137) 

K(Qi,P2)  =  -  ^  [a  -  c  cos  0Q1  -  bG  (P2/Qi)1  (138) 

The  integral  equations  then  become,  by  (105)  and  (106), 

4tcoi(Pi)  +  JSl  [oi(Qi)  -  ai(Pi)]  K(Pi,Qi)  dSi 
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(139) 


+  Js2  £°2(Q2)  K(Pi,Q2)  -  G2(p')  K(Q2,Pi)]  dSQ2  =  Ui 


9x 

dnpi 


4tCG2(P2)  +  [a2(Q2)  ■  G2(?2)}  K(P2/Q2)  dSQ2 

+  lSl  [oi(Ql)  K(P2,Qi)  -  oitpj,)  K(Qi,P2)]  dSqi  =  U2 (140) 

/  / 

in  which  0pi'  and  0p2'  at  Pj  and  P2  are  given  by  (117),  so  that  (139,  140)  are  in 

the  form  for  removing  both  the  singularities  of  the  kernels  and  the  peaks  of 
the  integrands. 

Because  of  the  axial  symmetry,  the  c's  are  not  functions  of  \| r.  We  may 
then  integrate  first  with  respect  to  xjtq,  and  since  the  integration  extends  over 
the  entire  circular  section,  we  may  take  yp  =  0.  For  bodies  of  revolution 
moving  in  the  direction  of  their  axes  of  symmetry,  it  is  known  that  the 
integration  over  the  polar  angle  y  in  the  plane  of  transverse  section  can  be 
expressed  in  terms  of  elliptic  integrals.  We  shall  take  advantage  of  this 
property  to  reduce  the  computing  effort  required  to  obtain  accurate  results  on 
a  3-D  problem,  and  will  treat  only  the  case  of  equal  spheres,  i.e.,  a  =  b. 

First,  then,  we  need  to  integrate  the  kernels  of  (139)  and  (140)  with 
respect  to  y.  For  the  cases  K(Pi,Qi)  and  K(P2,Q2),  consider  (with  i  =  1  or  2) 

1  -  G(Pj,Qi)  =  1  -  cos  0pi  cos  0Qj  -  sin  0pi  sin  0Qj  (2  cos2  ^  - 1) 

=  1  -  cos  (0pi  +  0Qi)  -  2  sin  0pi  sin  0Qi  sin2  v,  v  =  ^  (k  -  y) 

=  2  sin2  \  (0P  +  0q)  [1  -  £2(0p,0q)  sin2v],  £piQi  =  —  ^  §in  0Qi 

sin2  2  (0P  +  0Q) 


Then  by  (134), 
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(141) 


2n  |  !  k/2 

j  K(Pi,Qi)  dy  =  ~2  esc  ~  (0Pi  +  0Qi)  j  (1  -  ^PiQi  sin2v)‘1/2  dv 

0  a  z  Q 


2  (0Pi  +  0Pi)  K(^piQj) 


where  K(^p,Qi)  is  the  complete  elliptic  integral  of  the  first  kind. 

For  the  case  a  =  b,  the  other  four  kernels,  given  in  (135,  136, 
are  of  the  form 

a[l-G(0,0',v)]  -  c  cos  {^,} 

K(P'Q)  {2a2[l-G(0,0',v)]  -  2ac  (cos  0  +  cos  0')  +  c2}3/2 
where  0  =  0i  and  0'  =  jc  -  02,  and 


1-G(0,0',v)  =  1  +  cos  (0  -0')-2  sin  0  sin  0'  sin2v 
The  numerator  of  (142)  may  then  be  expressed  as 


,0' 


1  .00 


2^  (rpQ  +  2ac  cos  ^  }  -  c2)  =  ^  [rpQ  -  c(g  +  4a  sin2  ^  (0  ))1 


.  3 


and  since  the  denominator  is  rpQ,  we  obtain 


K(P,Q)  =  ^  ^  (g  +  4a  sin2 1 


PQ 


Also  one  can  show  that 


rPQ  =  D2[l  -  £2(0,0O  sin2  v] 


where 


D2  =  [g  +  2a(sin2  j  +  sin2  y)]2  +  a2(sin0  +  sin  0O2 


137,  138), 


(142) 


(143) 


(144) 


(145) 


(146a) 
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£2(8,0')  =  sin  8  sin  0' 


(146b) 


We  see  from  (144)  and  (145)  that  the  integral  of  the  first  term  of  the  kernel  is 
immediately  expressible  in  terms  of  an  elliptic  integral.  To  verify  that  the 
integral  of  the  second  term  is  also  so  expressible,  we  now  show  that 

71/2 

!  (1  -  $2)  (1 .  ^2  sin2V)-3/2  dv  =  E(5)  (147) 

0 

where  E(£)  is  the  complete  elliptic  integral  of  the  second  kind. 

The  integrand  of  (147)  may  be  written  as 


(l-£2  sin2v)_1/2  -  ^2cos2v(l-^2sin2v)‘3/2  (148) 


Integrating  the  second  term  by  parts,  we  obtain 
nl2  Till 

}  £2cos2v  (l-^2sin2v)‘3/2  dv  =  }  £2sin2v  (l-£2  sin2v)-!/2  dv 

which,  combined  with  the  integral  of  the  first  term  of  (148),  gives 

n/2 

E(^)  =  J  (1  -  ^sin^)1/2  dv 
0 


as  stated. 

The  integral  of  the  kernel  (144)  with  respect  to  \j/  when  P  and  Q  lie  on 
different  spheres  may  now  be  written  in  terms  of  elliptic  integrals.  We 
obtain,  by  (145), 


2n 


nil 


1  0' 

2c(g  +  4asin2  ^  {Q  }) 


!  K(P,Q)  dV  =  4  J  K(P,Q)  dv  =  s  KR«M»]  -  -  E©  04® 
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The  four  integrals  in  (139,  140)  can  now  be  evaluated  by  multiplying  the 
results  in  (141)  and  (149)  by  a2sin0Qi  or  a2sin0Q2  and  integrating  with  respect 

f 

to  0qi  or  0Q2  from  0  to  k.  These  become 


IPiQi  =  f  [Cfi(Qi)  -  C7j(Pi)]  sin  0q*  CSC  ^  (0Qi  +  6Pi)  K(^PiQi)  d0Qi;  1  =  (15°) 


n 


2a 


IP1Q2  *  J{[CT2(Q2)  -  02(P2)1  Dp^  K^P1Q2)  " 


acE(£,pi02) 

3  2 

DP1Q2  ^  ‘  ^P1Q2^ 


/  / 


[a2(Q2)  (g+4a  sin2  ^  0Q2 )  -  2o2(P2)  (g  +  4a  sin2  ^  0pi)])  sin  0Q2  d0Q2  (151) 


where,  by  (146), 


2  1  l  / 

DP1Q2  =  tg  +  2a  (sin2  2  0pi  +  sin2  2  0Q2^2  + 


and 


/  2 

a2 (sin  0pi  +  sin  0Q2)2  =  DQ2pi 


(152) 


2  sin0pjsin0Qj  K2  .  _ 

SpiQi  i  -  'kSiPi' 1  ~  L'z 

sin2^(0pi  +  0Qi) 

2  4a2  /  2 

^P1Q2  =  2  Sin  0pi  Sin  0Q2  =  ^Q2P1 

DP1Q2 


(153) 


(154) 


2  2 

Expressions  for  Ip2Qi,  Dp2Q!  and  ^P2Q!  can  be  obtained  from  (151,  152,  154)  by 
interchanging  subscripts  1  and  2. 
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K(£)  has  a  logarithmic  singularity  at  £  =  1  or  £'  =  (1  -  £2)1/2)  =  q.  When 
£'2  is  less  than  1/2,  K(£)  and  E(£)  can  be  computed  accurately  and  efficiently  by 
the  expansions 


K(£)  =  in  -  +  ^  4  (in  bn)  £'2n 


(155) 


E(£) 


4  1 


-  1  +  2  (in  ^  -  2 )  £'2  +  an  an+i  [in  ^  -  (2n+1)(2n+2)  '  bn]  ^ 


2n+2 


2n-l  1.3.5...(2n-l)  1  .1 

an  -  2n  an-i  -  2A  6  2n  ^  -  Vl  +  n(2n_D  -  ^  i(2i-l) 


(156) 


For  £2  less  than  1/2,  the  series  expansions 


(157) 


converge  rapidly. 

Since  the  integrands  F(0)  of  (150)  and  (151)  are  not  cyclic  in  0  <  0  <  k, 

f 

where  0  represents  0q  or  0q,  the  MAQF  cannot  be  used  directly  to  evaluate 

their  integrals.  We  could  generate  a  continuous  cyclic  function  by  defining 
F(0)  =  F(2jt-0)  for  k  <  Q  <  2n  and  then  repeating  this  combination  for 
successive  27c -intervals  to  the  left  and  right  of  the  0  to  2tc  range;  but  the 
derivatives  at  the  junction  points  0  =  0,  n,  2k  would  have  opposite  signs  since 
F’(0)  =  -  F(27C-0).  Hence  the  MAQF,  which  requires  a  smooth  function,  may 
not  be  applied  unless  the  slopes  are  zero  there. 

We  have  seen  that  the  nonuniform  MAQF  transformation  (107) 
introduces  the  factor  sin2  co/2  in  the  integrand,  which  gives  zero  derivatives 
with  respect  to  co  at  co  =  0  and  2k,  but  not  at  co  =  k.  By  replacing  0  and  to  in  (107) 
by  20  and  2co,  the  transformation  and  its  differential  become 

1 

0  =  to  -  2  sin  2co,  d0  =  2sin2  tod  to  (158) 


and  then 
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(159) 


k  n 

j  F(0)  de  =  2  J  F(6)  sin2  todco 
0  0 

Since  sin  0  is  a  factor  of  F(0),  the  zeros  of  the  integrand  of  (159)  at  co  =  0,  k,  2k 
are  seen  to  be  of  third  order.  Hence  not  only  the  slopes  but  also  the 
curvatures  of  the  extended  integrands  are  zero  at  the  junction  points  as  a 
consequence  of  the  transformation  (158).  If  a  higher  concentration  of  points 
near  0  =  0  than  that  given  by  (158)  is  required,  the  transformations  (107a)  and 
(107b)  may  be  used,  again  with  0  and  co  replaced  by  20  and  2co. 

Results  for  ki,  ki2,  kj,  ki2  are  given  in  Table  12,  where  calculated  values 
are  compared  with  the  exact  ones.  The  derivatives  were  obtained  by  the  same 
procedure  as  for  the  circles,  i.e.,  numerical  differentiation  with  respect  to  the 
parameter  y  by  a  five  point  central-difference  formula,  followed  by  the  exact 
derivative  of  y  with  respect  to  g,  as  is  shown  in  (109).  The  results  for  ki  and 
ki2  practically  agree  to  five  decimals  from  g/b  =  0.01  to  0.04,  and  to  six 
decimals  at  larger  gaps;  ki  and  ki2  agree  to  three  decimals  up  to  g/b  =  0.03  to 
four  decimals  up  to  0.07,  and  to  five  decimals  at  larger  gaps. 

Results  for  kj  and  k12  are  compared  with  those  of  Ref.  [2]  in  Table  13. 
The  deviation  from  the  'exact'  values  is  seen  to  be  much  less  with  the  present 
procedures,  which  were  accurate  enough  for  application  of  numerical 
differentiation  to  obtain  the  added-mass  derivatives  given  in  Table  12.  Only 
exact  values  of  these  derivatives  are  presented  in  Ref.  [2],  indicating  that  their 
procedure  of  solving  the  pair  of  differentiated  integral  equations  could  not 
give  adequate  accuracy  at  small  gaps,  probably  because  the  peaks  of  their 
differentiated  kernels  are  of  the  order  of  1/g3. 

SUMMARY  AND  CONCLUSIONS 

The  motivation  for  undertaking  the  present  work  was  to  evaluate  the 
accuracy  of  the  interaction  forces  computed  by  the  integral-equation  approach. 
For  this  purpose,  the  classical  theory  of  successive  images,  for  central  impact 
of  pairs  of  circles  or  spheres,  is  available  for  obtaining  highly  accurate 
numerical  solutions,  against  which  results  from  integral  equations  could  be 
compared.  In  the  course  of  applying  these  approaches,  new  and  significant 
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advances  in  both  areas  were  made  and  have  been  described.  These  will  now 
be  summarized. 

On  the  method  of  successive  images,  the  following  has  been 
accomplished: 

1.  A  new  and  more  rational  derivation  of  Herman’s  formulas  for  the 
doublet  strengths  than  that  given  in  his  paper  [10]  is  presented. 

2.  A  parametric  form  of  these  formulas  has  been  applied  to  derive  new 
truncation-correction  formulas,  with  which  the  infinite  series  for  the  added 
masses  can  be  summed  to  a  desired  accuracy  with  many  fewer  terms. 

3.  A  new  asymptotic  formula  equating  the  mathematical  parameter  £ 
of  the  aforementioned  parametric  form  to  a  polynomial  in  the 
nondimensional  physical  parameter  y  (which  is  proportional  to  the  square- 
root  of  the  gap  g  between  the  bodies)  is  derived,  and  shown  to  be  highly 
accurate  at  small  gaps,  and  for  equal  circular  cylinders  or  spheres,  to  be 
usefully  accurate  up  to  g/b  =  1. 

4.  Combining  the  truncation  and  asymptotic  formulas  yielded  the 
results  that  the  series  for  the  added  masses,  and  those  for  their  derivatives 
with  respect  to  the  parameter  £,  both  converge  at  g  =  0,  and  asymptotic  Taylor 
expansions  about  g  =  0,  which  display  the  added  masses  and  their  derivatives 

y  at  g  =  0,  are  given.  These  properties  led  to  the  new  result  that  the  derivatives 
of  the  added  masses  with  respect  to  g,  for  small  values  of  y,  vary  as  gr1/2  for 
circular  cylinders  and  as  -in  g  for  spheres.  In  both  cases,  this  implies  that  the 
repulsive  forces  approach  infinity,  verifying  the  well-known  irrotational-flow 
paradox  that  the  bodies  would  never  meet. 

5.  The  need  for  considering  the  property  of  uniform  convergence  of 
the  series  for  the  added  masses  and  their  derivatives  is  discussed  in  the  text. 
It  is  proved  that  these  series  are  uniformly  convergent  in  the  closed  region  0 
<  g  <  °°>  indicating  that  the  added  masses  are  continuous  functions  of  g  at  g  = 
0;  and  that  the  series  for  their  derivatives  with  respect  to  £  are  uniformly 
convergent  in  the  open  region  0  <  g  <  <»,  i.e.,  although  this  derivative  series 
converges  at  g  =  0,  the  convergence  is  not  uniform.  This  implies  that  the 


60 


y  derivatives  of  the  added  masses  may  be  discontinuous  at  £  =  0,  as  is  shown  in 
^  the  text  for  a  pair  of  circular  cylinders. 

On  the  method  of  integral  equations,  three  new  procedures  were 
developed  to  obtain  more  accurate  solutions.  The  first  of  these  was  required 
at  small  gaps  to  eliminate  the  peaks  of  the  four  kernels  of  the  two  integral 
equations.  It  was  shown  that  the  transposes  of  these  kernels  eliminate  not 
only  their  singularities  but  also  their  peaks,  which  are  of  the  order  of  1  /g  for 
two-dimensional  and  of  1/g2  for  three-dimensional  bodies  of  general  shape. 
This  left  a  residue  of  much  smaller  peaks  due  to  the  rapid  variation  of  the 
source  distribution  in  the  gap  region  when  the  gap  is  small.  These  were 
treated  by  applying  a  quadrature  formula,  described  in  the  next  paragraph, 
which  concentrates  many  points  in  the  small  neighborhood  of  the  smaller 
peaks. 

In  the  second  procedure,  the  'most  accurate  quadrature  formula,' 
which  requires  a  smooth,  cyclic  integrand  and  uniform  intervals,  was 
modified  by  changing  the  variable  of  integration,  so  that  the  integrand 
remains  cyclic,  uniform  intervals  are  taken  in  the  new  variable,  and  many 
points  of  the  original  variable  are  concentrated  in  the  desired  region.  A 
sequence  of  such  transformations,  of  successively  increasing  point 
concentrations  in  a  small  region,  is  presented.  With  a  slight  modification, 
this  nonuniform  MAQF  also  served  to  evaluate  accurately  the  noncyclic 
integrals  of  the  two-sphere  problem. 

The  third  procedure  was  developed  in  order  to  improve  the  accuracy  of 
the  added-mass  derivatives  obtained  by  numerical  differentiation  of  the 
added-mass  data.  As  is  described  in  the  text,  this  procedure,  which  is  suitable 
for  general  shapes,  requires  accurate  solutions  of  the  integral  equations,  and  is 
suggested  by  the  asymptotic  formulas  for  the  added  masses  at  small  gaps. 

By  means  of  these  new  procedures,  the  minimum  gap  at  which 
accurate  results  could  be  obtained  was  reduced  to  about  one-tenth  of  those 
reported  in  the  original  work  on  this  problem  for  the  Mobil  Research  and 
Development  Corporation.  Results  for  circle  pairs  of  diameter  ratios  1,4,16 
and  «>,  two  ellipse-circle  pairs  and  a  pair  of  equal  spheres  are  presented.  The 
treatment  of  the  infinite-diameter  ratio  (a  circle  approaching  a  wall),  which  is 
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also  applied  to  a  sphere  and  a  wall,  yields  results  for  the  added  masses  and 
forces  on  bodies  of  arbitrary  shape.  All  the  other  applications  were  to  two- 
dimensional  and  axisymmetric  forms,  for  which  the  integral-equation 
approach  yielded  highly  accurate  solutions.  An  IBM  RISC/6000  was  used  to 
obtain  the  numerical  results. 
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Table  1.  Example  showing  truncation  correction  by  trapezoidal  integral 


N 

R(N) 

N  1 

L,  ^  +  R(N) 

10 

0.095238 

1.645006 

100 

0.009950 

1.644934 

1000 

0.001000 

1.644934 

Table  2.  Calculation  of  -k\2  with  truncation  correction  R12  for  a=b. 


g/a  = 

0.0001 

Zl 

g/a 

N 

R12 

-ki2** 

Rl2 

-ki2 

Rl2 

-ki2 

10 

.0383183 

.8125634 

.0228071 

.7914781 

.0030591 

.7282811 

20 

.0157424 

.8125323 

.0051146 

.7914469 

.0000553 

.7282603 

30 

.0083784 

.8125290 

.0013646 

.7914436 

.0000010 

.7282599 

40 

.0049347 

.8125281 

.0003793 

.7914429 

0 

• 

50 

.0030590 

.8125278 

.0001066 

.7914427 

. 

• 

70 

.0012677 

.8125276 

.0000085 

. 

. 

• 

100 

.0003656 

. 

.0000002 

. 

• 

• 

200 

.0000066 

. 

. 

. 

. 

• 

300 

.0000001 

. 

• 

. 

• 

OO 

0 

.8125275 

0 

.7914426 

0 

.7282599 

*  -kl2=£>£^82n-l  +R12 


Table  3.  Values  of  ki  by  various  methods  at  g/b  =  0.01  for  a  =  b 


ki 

Inteerand  peaks 

Nonuniform 

N  =  40 

N  =  100 

removed? 

MAOF 

1. 

Yes 

Yes 

1.375338 

1.375338 

2. 

No.  onlv  kernel 

Yes 

1.0177 

1.3715 

peaks  removed* 

3. 

Yes 

uniform 

1.3790 

1.3785 

4. 

No 

Yes 

6.8110 

1.3748 

5. 

No 

uniform 

diverged 

diverged 

*Used  ai(0)  and  02(11)  instead  of  oi(0ip')  and  O2(02P') 
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Table.  4.  Comparison  of  'exact'  added  masses  and  their  derivatives  for  equal 
circles  with  results  calculated  using  integral  equations. 

k^calc)  k^(calc)  k^(exact)  k  ^  2<ca*c)  ki2(ca*c)  ki2(exact)  (-'F1"'<-'F2  Cp2»-C| 

1.375338  -4.22054  -4.22048  -0.728260  4.43439  4.43434  -6.6296  7.3014 
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Table  5.  Comparison  with  Ref  [2]  of  values  of  ki  and  k\2  for  equal  circles. 


ki  _ -ki2 


g/b 

exact 

present* 
int.  eq. 

Ref.  [2]** 
int.  eq. 

exact 

present 
int.  eq. 

Ref.  [2] 
int.  eq. 

0.01 

1.375338 

1.375338 

1.3713 

0.728260 

0.728260 

0.7214 

0.02 

1.341620 

1.341620 

1.3400 

0.692409 

0.692409 

0.6901 

0.03 

1.317390 

1.317390 

1.3169 

0.666057 

0.666057 

0.6653 

0.04 

1.298030 

1.298030 

1.2979 

0.644586 

0.644586 

0.6443 

0.05 

1.281755 

1.281755 

1.2817 

0.626212 

0.626212 

0.6261 

0.06 

1.267655 

1.267655 

1.2676 

0.610023 

0.610023 

0.6100 

0.07 

1.255188 

1.255188 

1.2552 

0.595480 

0.595480 

0.5955 

0.08 

1.244005 

1.244005 

1.2440 

0.582231 

0.582231 

0.5822 

0.09 

1.233862 

1.233862 

1.2339 

0.570034 

0.570034 

0.5700 

0.10 

1.224583 

1.224583 

1.2246 

0.558711 

0.558711 

0.5587 

^Calculated  with  a  40-point  nonuniform  MAQF,  and  all  integrand  peaks  removed. 


'"'"Calculated  with  three  Gauss  quadrature  formulas,  of  order  20  (see  condition  2  of  Table  3). 


Table  6.  Comparison  with  Ref  [2]  of  values  of  ki  and  ki2  for  a  =  b  =  1. 
_  -ki  _ki2 _ 


g/b 

exact 

present* 
int.  eq. 

exact 

present* 
int.  eq. 

0.01 

4.220480 

4.22054 

4.10983 

4.434339 

4.43439 

4.32371 

0.02 

2.773348 

2.77335 

2.77325 

2.986088 

2.98609 

2.98597 

0.03 

2.137694 

2.13769 

2.13769 

2.349315 

2.34931 

2.34931 

0.04 

1.761898 

1.76190 

1.76190 

1.972399 

1.97240 

1.97240 

0.05 

1.507534 

1.50753 

1.50753 

1.716914 

1.71691 

1.71691 

0.06 

1.321292 

1.32129 

1.32129 

1.529551 

1.52955 

1.52955 

0.07 

1.177713 

1.17771 

1.17771 

1.384850 

1.38485 

1.38485 

0.08 

1.062909 

1.06291 

1.06291 

1.268923 

1.26892 

1.26892 

0.09 

0.968580 

0.96858 

0.96868 

1.173472 

1.17347 

1.17347 

0.10 

0.889422 

0.88942 

0.88942 

1.093190 

1.09319 

1.09319 

■"Calculated  from  a  five-point  central-difference  formula,  for  which  values  of  ki  and 
ki2  at  g/b  =  0.005  and  0.007  were  required  and  computed. 
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Table  7.  Results  for  unequal  circles; 
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Table  8.  Results  for  unequal  circles;  a  =  16b. 
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Table  9.  Results  for  a  circle  and  a  wall;  a  =  «>. 
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Table  10.  Results  for  an  ellipse  and  a  circle;  a,  =  2b,  ^  = 
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Table  11.  Results  for  an  ellipse  and  a  circle;  aj  =  b,  a2  =  2b. 
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Table.  12.  Comparison  of  'exact1  and  calculated  values  of  added  masses 
and  their  derivatives  for  equal  spheres. 
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Table  13.  Comparison  with  Ref.  [2]  of  values  of  ki  and  ki2  for  equal 

spheres. 


ki 

-ki2 

g/b 

'exact' 

calculated 

Ref.  [2]  calc. 

'exact' 

calculated 

Ref.  [2]  calc. 

0.01 

0.570162 

0.570146 

— 

0.216292 

0.216276 

— 

0.02 

0.565138 

0.565129 

0.5614 

0.209728 

0.209719 

0.2063 

0.03 

0.561022 

0.561018 

0.5591 

0.204085 

0.204080 

0.2025 

0.04 

0.557477 

0.557474 

0.5562 

0.199028 

0.199025 

0.1981 

0.05 

0.554343 

0.554341 

0.5540 

0.194396 

0.194394 

0.1943 

0.06 

0.551527 

0.551526 

0.5512 

0.190098 

0.190097 

0.1900 

0.07 

0.548968 

0.548967 

0.5486 

0.186071 

0.186070 

0.1860 

0.08 

0.546625 

0.546625 

0.5464 

0.182273 

0.182273 

0.1823 

0.09 

0.544464 

0.544464 

- 

0.178672 

0.178672 

- 

0.10 

0.542461 

0.542461 

0.5421 

0.175244 

0.175244 

0.1752 
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Figure  1 .  Definition  sketch  of  two  circles  or  spheres  moving  along  then- 
line  of  centers  for  method  of  successive  images. 
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Figure  3.  Sketch  of  two  cylinders  for  method  of  integral  equations. 


Figure  4.  Dlustration  of  procedure  for  eliminating  peaks  of  integrand. 
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moving  cylinder  a  moving  cylinder  D  moving  cylinder 

stationary  cylinder  ♦  stationary  cylinder  ♦  stationary  cylinder 
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Figure  5.  Source  distributions  on  equal  circles  for  Uj  =  1  and  U. 


Figure  6.  Integrands  of  moving  circle  for  various  fixed  points  on 
stationary  one  for  a  =  b  and  g/b  =  0.01. 
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